CVC3
rational-native.cpp
Go to the documentation of this file.
1 /*****************************************************************************/
2 /*!
3  * \file rational-native.cpp
4  *
5  * \brief Implementation of class Rational using native (bounded
6  * precision) computer arithmetic
7  *
8  * Author: Sergey Berezin
9  *
10  * Created: Mon Jul 28 12:18:03 2003
11  *
12  * <hr>
13  * License to use, copy, modify, sell and/or distribute this software
14  * and its documentation for any purpose is hereby granted without
15  * royalty, subject to the terms and conditions defined in the \ref
16  * LICENSE file provided with this distribution.
17  *
18  * <hr>
19  *
20  */
21 /*****************************************************************************/
22 
23 #ifdef RATIONAL_NATIVE
24 
25 #include "compat_hash_set.h"
26 #include "rational.h"
27 // For atol() (ASCII to long)
28 #include <stdlib.h>
29 #include <limits.h>
30 #include <errno.h>
31 #include <sstream>
32 #include <math.h>
33 
34 #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."
35 
36 namespace CVC3 {
37 
38  using namespace std;
39 
40  //! Add two integers and check for overflows
41  static long int plus(long int x, long int y) {
42  long int res = x+y;
43  FatalAssert(((x > 0) != (y > 0)) || ((x > 0) == (res > 0)),
44  "plus(x,y): arithmetic overflow" OVERFLOW_MSG);
45  return res;
46  }
47 
48  //! Add two integers and check for overflows
49  static unsigned long uplus(unsigned long x, unsigned long y) {
50  unsigned long res = x+y;
51  FatalAssert(res >= x && res >= y,
52  "uplus(x,y): arithmetic overflow" OVERFLOW_MSG);
53  return res;
54  }
55 
56  //! Subtract two unsigned integers and check for overflows
57  static unsigned long unsigned_minus(unsigned long x, unsigned long y) {
58  unsigned long res = x-y;
59  FatalAssert(res <= x,
60  "unsigned_minus(x,y): arithmetic overflow" OVERFLOW_MSG);
61  return res;
62  }
63 
64  //! Multiply two unsigned integers and check for overflows
65  static unsigned long umult(unsigned long x, unsigned long y) {
66  unsigned long res = x*y;
67  FatalAssert(res == 0 || res/x == y,
68  "umult(x,y): arithmetic overflow" OVERFLOW_MSG);
69  return res;
70  }
71 
72  //! Shift two unsigned integers and check for overflow
73  static unsigned long ushift(unsigned long x, unsigned y) {
74  FatalAssert(y < (unsigned)numeric_limits<unsigned long>::digits,
75  "ushift(x,y): arithmetic overflow" OVERFLOW_MSG);
76  unsigned long res = (x << y);
77  FatalAssert((res >> y) == x,
78  "ushift(x,y): arithmetic overflow" OVERFLOW_MSG);
79  return res;
80  }
81 
82  //! Compute GCD using Euclid's algorithm (from Aaron Stump's code)
83  static unsigned long ugcd(unsigned long n1, unsigned long n2) {
84  DebugAssert(n1!=0 && n2!=0,
85  "gcd("+int2string(n1)+", "+int2string(n2)+"): bad args");
86  if (n1 < n2) {
87  long int tmp = n1;
88  n1 = n2;
89  n2 = tmp;
90  }
91 
92  // at this point, n1 >= n2
93  long int r = n1 % n2;
94  if (!r)
95  return n2;
96 
97  return ugcd(n2, r);
98 }
99 
100  //! Unary minus which checks for overflows
101  static long int uminus(long int x) {
102  FatalAssert(x == 0 || x != -x, "uminus(x): arithmetic overflow"
103  OVERFLOW_MSG);
104  return -x;
105  }
106 
107  //! Multiply two integers and check for overflows
108  static long int mult(long int x, long int y) {
109  long int res = x*y;
110  FatalAssert(x==0 || res/x == y, "mult(x,y): arithmetic overflow"
111  OVERFLOW_MSG);
112  return res;
113  }
114 
115  //! Compute GCD using Euclid's algorithm (from Aaron Stump's code)
116  static long int gcd(long int n1, long int n2) {
117  DebugAssert(n1!=0 && n2!=0,
118  "gcd("+int2string(n1)+", "+int2string(n2)+"): bad args");
119  // First, make the arguments positive
120  if(n1 < 0) n1 = uminus(n1);
121  if(n2 < 0) n2 = uminus(n2);
122 
123  if (n1 < n2) {
124  long int tmp = n1;
125  n1 = n2;
126  n2 = tmp;
127  }
128 
129  // at this point, n1 >= n2
130  long int r = n1 % n2;
131  if (!r)
132  return n2;
133 
134  return gcd(n2, r);
135 }
136 
137  //! Compute LCM
138  static long int lcm(long int n1, long int n2) {
139  long int g = gcd(n1,n2);
140  return mult(n1/g, n2);
141  }
142 
143  static long int ulcm(unsigned long n1, unsigned long n2) {
144  long int g = ugcd(n1,n2);
145  return umult(n1/g, n2);
146  }
147 
148  // Implementation of the forward-declared internal class
149  class Rational::Impl {
150  long int d_num; //!< Numerator
151  long int d_denom; //!< Denominator
152  //! Make the rational number canonical
153  void canonicalize();
154  public:
155  //! Default Constructor
156  Impl(): d_num(0), d_denom(1) { }
157  //! Copy constructor
158  Impl(const Impl &x) : d_num(x.d_num), d_denom(x.d_denom) { }
159  //! Constructor from unsigned long
160  Impl(unsigned long n) :d_num(n), d_denom(1) {
161  FatalAssert(d_num >= 0, "Rational::Impl(unsigned long) - arithmetic overflow");
162  }
163  //! Constructor from a pair of integers
164  Impl(long int n, long int d): d_num(n), d_denom(d) { canonicalize(); }
165  //! Constructor from a string
166  Impl(const string &n, int base);
167  //! Constructor from a pair of strings
168  Impl(const string &n, const string& d, int base);
169  // Destructor
170  virtual ~Impl() { }
171  //! Get numerator
172  long int getNum() const { return d_num; }
173  //! Get denominator
174  long int getDen() const { return d_denom; }
175 
176  //! Unary minus
177  Impl operator-() const;
178 
179  //! Equality
180  friend bool operator==(const Impl& x, const Impl& y) {
181  return (x.d_num == y.d_num && x.d_denom == y.d_denom);
182  }
183  //! Dis-equality
184  friend bool operator!=(const Impl& x, const Impl& y) {
185  return (x.d_num != y.d_num || x.d_denom != y.d_denom);
186  }
187  /*!
188  * Compare x=n1/d1 and y=n2/d2 as n1*f2 < n2*f1, where f1=d1/f,
189  * f2=d2/f, and f=lcm(d1,d2)
190  */
191  friend bool operator<(const Impl& x, const Impl& y) {
192  Impl diff(x-y);
193  return diff.d_num < 0;
194  }
195 
196  friend bool operator<=(const Impl& x, const Impl& y) {
197  return (x == y || x < y);
198  }
199 
200  friend bool operator>(const Impl& x, const Impl& y) {
201  Impl diff(x-y);
202  return diff.d_num > 0;
203  }
204 
205  friend bool operator>=(const Impl& x, const Impl& y) {
206  return (x == y || x > y);
207  }
208 
209  /*! Addition of x=n1/d1 and y=n2/d2: n1*g2 + n2*g1, where g1=d1/g,
210  * g2=d2/g, and g=lcm(d1,d2)
211  */
212  friend Impl operator+(const Impl& x, const Impl& y) {
213  long int d1(x.getDen()), d2(y.getDen());
214  long int f(lcm(d1,d2)), f1(f/d1), f2(f/d2);
215  long int n = plus(mult(x.getNum(), f1), mult(y.getNum(), f2));
216  return Impl(n, f);
217  }
218 
219  friend Impl operator-(const Impl& x, const Impl& y) {
220  TRACE("rational", "operator-(", x, ", "+y.toString()+")");
221  long int d1(x.getDen()), d2(y.getDen());
222  long int f(lcm(d1,d2)), f1(f/d1), f2(f/d2);
223  long int n = plus(mult(x.getNum(), f1), uminus(mult(y.getNum(), f2)));
224  Impl res(n, f);
225  TRACE("rational", " => ", res, "");
226  return res;
227  }
228 
229  /*!
230  * Multiplication of x=n1/d1, y=n2/d2:
231  * (n1/g1)*(n2/g2)/(d1/g2)*(d2/g1), where g1=gcd(n1,d2) and
232  * g2=gcd(n2,d1)
233  */
234  friend Impl operator*(const Impl& x, const Impl& y) {
235  long int n1(x.getNum()), d1(x.getDen()), n2(y.getNum()), d2(y.getDen());
236  long int g1(n1? gcd(n1,d2) : 1), g2(n2? gcd(n2,d1) : 1);
237  long int n(mult(n1/g1, n2/g2)), d(mult(d1/g2, d2/g1));
238  return Impl(n,d);
239  }
240 
241  /*!
242  * Division of x=n1/d1, y=n2/d2:
243  * (n1/g1)*(d2/g2)/(d1/g2)*(n2/g1), where g1=gcd(n1,n2) and
244  * g2=gcd(d1,d2)
245  */
246  friend Impl operator/(const Impl& x, const Impl& y) {
247  long int n1(x.getNum()), d1(x.getDen()), n2(y.getNum()), d2(y.getDen());
248  DebugAssert(n2 != 0, "Impl::operator/: divisor is 0");
249  long int g1(n1? gcd(n1,n2) : 1), g2(gcd(d1,d2));
250  long int n(n1? mult(n1/g1, d2/g2) : 0), d(n1? mult(d1/g2, n2/g1) : 1);
251  return Impl(n,d);
252  }
253 
254  friend Impl operator%(const Impl& x, const Impl& y) {
255  DebugAssert(x.getDen() == 1 && y.getDen() == 1,
256  "Impl % Impl: x and y must be integers");
257  return Impl(x.getNum() % y.getNum(), 1);
258  }
259 
260  //! Print to string
261  string toString(int base = 10) const {
262  ostringstream ss;
263  if (d_num == 0) ss << "0";
264  else if (base == 10) {
265  ss << d_num;
266  if (d_denom != 1)
267  ss << "/" << d_denom;
268  }
269  else {
270  vector<int> vec;
271  long num = d_num;
272  while (num) {
273  vec.push_back(num % base);
274  num = num / base;
275  }
276  while (!vec.empty()) {
277  if (base > 10 && vec.back() > 10) {
278  ss << (char)('A' + (vec.back()-10));
279  }
280  else ss << vec.back();
281  vec.pop_back();
282  }
283  if(d_denom != 1) {
284  ss << "/";
285  if (d_denom == 0) ss << "0";
286  else {
287  num = d_denom;
288  while (num) {
289  vec.push_back(num % base);
290  num = num / base;
291  }
292  while (!vec.empty()) {
293  if (base > 10 && vec.back() > 10) {
294  ss << (char)('A' + (vec.back()-10));
295  }
296  else ss << vec.back();
297  vec.pop_back();
298  }
299  }
300  }
301  }
302  return(ss.str());
303  }
304 
305  //! Printing to ostream
306  friend ostream& operator<<(ostream& os, const Rational::Impl& n) {
307  return os << n.toString();
308  }
309 
310  };
311 
312  // Make the rational number canonical
313  void Rational::Impl::canonicalize() {
314  DebugAssert(d_denom != 0,
315  "Rational::Impl::canonicalize: bad denominator: "
316  +int2string(d_denom));
317  if(d_num == 0) {
318  d_denom = 1;
319  } else {
320  if(d_denom < 0) {
321  d_num = uminus(d_num);
322  d_denom = uminus(d_denom);
323  }
324  long int d = gcd(d_num, d_denom);
325  if(d != 1) {
326  d_num /= d;
327  d_denom /= d;
328  }
329  }
330  }
331 
332  // Constructor from a string
333  Rational::Impl::Impl(const string &n, int base) {
334  size_t i, iend;
335  for(i=0,iend=n.size(); i<iend && n[i] != '/'; ++i);
336  if(i<iend)
337  // Found slash at i
338  *this = Impl(n.substr(0,i-1), n.substr(i+1,iend-1), base);
339  else
340  *this = Impl(n, "1", base);
341  canonicalize();
342  }
343  // Constructor from a pair of strings
344  Rational::Impl::Impl(const string &n, const string& d, int base) {
345  d_num = strtol(n.c_str(), NULL, base);
346  FatalAssert(d_num != LONG_MIN && d_num != LONG_MAX,
347  "Rational::Impl(string, string): arithmetic overflow:"
348  "n = "+n+", d="+d+", base="+int2string(base)
349  +OVERFLOW_MSG);
350  d_denom = strtol(d.c_str(), NULL, base);
351  FatalAssert(d_denom != LONG_MIN && d_denom != LONG_MAX,
352  "Rational::Impl(string, string): arithmetic overflow:"
353  "n = "+n+", d="+d+", base="+int2string(base)
354  +OVERFLOW_MSG);
355  canonicalize();
356  }
357 
358  Rational::Impl Rational::Impl::operator-() const {
359  return Impl(uminus(d_num), d_denom);
360  }
361 
362 
363  //! Default constructor
364  Rational::Rational() : d_n(new Impl) {
365 #ifdef _DEBUG_RATIONAL_
366  int &num_created = getCreated();
367  num_created++;
368  printStats();
369 #endif
370  }
371  //! Copy constructor
372  Rational::Rational(const Rational &n) : d_n(new Impl(*n.d_n)) {
373 #ifdef _DEBUG_RATIONAL_
374  int &num_created = getCreated();
375  num_created++;
376  printStats();
377 #endif
378  }
379 
380  Rational::Rational(const Unsigned &n) : d_n(new Impl(n.getUnsigned())) {
381 #ifdef _DEBUG_RATIONAL_
382  int &num_created = getCreated();
383  num_created++;
384  printStats();
385 #endif
386  }
387 
388  //! Private constructor
389  Rational::Rational(const Impl& t): d_n(new Impl(t)) {
390 #ifdef _DEBUG_RATIONAL_
391  int &num_created = getCreated();
392  num_created++;
393  printStats();
394 #endif
395  }
396  Rational::Rational(int n, int d): d_n(new Impl(n, d)) {
397 #ifdef _DEBUG_RATIONAL_
398  int &num_created = getCreated();
399  num_created++;
400  printStats();
401 #endif
402  }
403  // Constructors from strings
404  Rational::Rational(const char* n, int base)
405  : d_n(new Impl(string(n), base)) {
406 #ifdef _DEBUG_RATIONAL_
407  int &num_created = getCreated();
408  num_created++;
409  printStats();
410 #endif
411  }
412  Rational::Rational(const string& n, int base)
413  : d_n(new Impl(n, base)) {
414 #ifdef _DEBUG_RATIONAL_
415  int &num_created = getCreated();
416  num_created++;
417  printStats();
418 #endif
419  }
420  Rational::Rational(const char* n, const char* d, int base)
421  : d_n(new Impl(string(n), string(d), base)) {
422 #ifdef _DEBUG_RATIONAL_
423  int &num_created = getCreated();
424  num_created++;
425  printStats();
426 #endif
427  }
428  Rational::Rational(const string& n, const string& d, int base)
429  : d_n(new Impl(n, d, base)) {
430 #ifdef _DEBUG_RATIONAL_
431  int &num_created = getCreated();
432  num_created++;
433  printStats();
434 #endif
435  }
436  // Destructor
437  Rational::~Rational() {
438  delete d_n;
439 #ifdef _DEBUG_RATIONAL_
440  int &num_deleted = getDeleted();
441  num_deleted++;
442  printStats();
443 #endif
444  }
445 
446  // Get components
447  Rational Rational::getNumerator() const
448  { return Rational(Impl(d_n->getNum(), 1)); }
449  Rational Rational::getDenominator() const
450  { return Rational(Impl(d_n->getDen(), 1)); }
451 
452  bool Rational::isInteger() const { return 1 == d_n->getDen(); }
453 
454  // Assignment
455  Rational& Rational::operator=(const Rational& n) {
456  if(this == &n) return *this;
457  delete d_n;
458  d_n = new Impl(*n.d_n);
459  return *this;
460  }
461 
462  ostream &operator<<(ostream &os, const Rational &n) {
463  return(os << n.toString());
464  }
465 
466 
467  // Check that argument is an int and print an error message otherwise
468 
469  static void checkInt(const Rational& n, const string& funName) {
470  DebugAssert(n.isInteger(),
471  ("CVC3::Rational::" + funName
472  + ": argument is not an integer: " + n.toString()).c_str());
473  }
474 
475  /* Computes gcd and lcm on *integer* values. Result is always a
476  positive integer. In this implementation, it is guaranteed by
477  GMP. */
478 
479  Rational gcd(const Rational &x, const Rational &y) {
480  checkInt(x, "gcd(*x*,y)");
481  checkInt(y, "gcd(x,*y*)");
482  return Rational(Rational::Impl(gcd(x.d_n->getNum(), y.d_n->getNum()), 1));
483  }
484 
485  Rational gcd(const vector<Rational> &v) {
486  long int g(1);
487  if(v.size() > 0) {
488  checkInt(v[0], "gcd(vector<Rational>[0])");
489  g = v[0].d_n->getNum();
490  }
491  for(size_t i=1; i<v.size(); i++) {
492  checkInt(v[i], "gcd(vector<Rational>)");
493  if(g == 0) g = v[i].d_n->getNum();
494  else if(v[i].d_n->getNum() != 0)
495  g = gcd(g, v[i].d_n->getNum());
496  }
497  return Rational(Rational::Impl(g,1));
498  }
499 
500  Rational lcm(const Rational &x, const Rational &y) {
501  long int g;
502  checkInt(x, "lcm(*x*,y)");
503  checkInt(y, "lcm(x,*y*)");
504  g = lcm(x.d_n->getNum(), y.d_n->getNum());
505  return Rational(Rational::Impl(g, 1));
506  }
507 
508  Rational lcm(const vector<Rational> &v) {
509  long int g(1);
510  for(size_t i=0; i<v.size(); i++) {
511  checkInt(v[i], "lcm(vector<Rational>)");
512  if(v[i].d_n->getNum() != 0)
513  g = lcm(g, v[i].d_n->getNum());
514  }
515  return Rational(Rational::Impl(g,1));
516  }
517 
518  Rational abs(const Rational &x) {
519  long int n(x.d_n->getNum());
520  if(n>=0) return x;
521  return Rational(Rational::Impl(-n, x.d_n->getDen()));
522  }
523 
524  Rational floor(const Rational &x) {
525  if(x.d_n->getDen() == 1) return x;
526  long int n = x.d_n->getNum();
527  long int nAbs = (n<0)? uminus(n) : n;
528  long int q = nAbs / x.d_n->getDen();
529  if(n < 0) q = plus(uminus(q), -1);
530  return Rational(Rational::Impl(q,1));
531  }
532 
533  Rational ceil(const Rational &x) {
534  if(x.d_n->getDen() == 1) return x;
535  long int n = x.d_n->getNum();
536  long int nAbs = (n<0)? -n : n;
537  long int q = nAbs / x.d_n->getDen();
538  if(n > 0) q = plus(q, 1);
539  else q = uminus(q);
540  return Rational(Rational::Impl(q,1));
541  }
542 
543  Rational mod(const Rational &x, const Rational &y) {
544  checkInt(x, "mod(*x*,y)");
545  checkInt(y, "mod(x,*y*)");
546  long int r = x.d_n->getNum() % y.d_n->getNum();
547  return(Rational(Rational::Impl(r,1)));
548  }
549 
550  Rational intRoot(const Rational& base, unsigned long int n) {
551  checkInt(base, "intRoot(*x*,y)");
552  checkInt(n, "intRoot(x,*y*)");
553  double b = base.d_n->getNum();
554  double root = n;
555  root = 1/root;
556  b = ::pow(b, root);
557  long res = (long) ::floor(b);
558  if (::pow((long double)res, (int)n) == base.d_n->getNum()) {
559  return Rational(Rational::Impl(res,1));
560  }
561  return Rational(Rational::Impl((long)0,(long)1));
562  }
563 
564  string Rational::toString(int base) const {
565  return(d_n->toString(base));
566  }
567 
568  size_t Rational::hash() const {
570  return h(toString().c_str());
571  }
572 
573  void Rational::print() const {
574  cout << (*this) << endl;
575  }
576 
577  // Unary minus
578  Rational Rational::operator-() const {
579  return Rational(Rational::Impl(-(d_n->getNum()), d_n->getDen()));
580  }
581 
582  Rational &Rational::operator+=(const Rational &n2) {
583  *this = (*this) + n2;
584  return *this;
585  }
586 
587  Rational &Rational::operator-=(const Rational &n2) {
588  *this = (*this) - n2;
589  return *this;
590  }
591 
592  Rational &Rational::operator*=(const Rational &n2) {
593  *this = (*this) * n2;
594  return *this;
595  }
596 
597  Rational &Rational::operator/=(const Rational &n2) {
598  *this = (*this) / n2;
599  return *this;
600  }
601 
602  int Rational::getInt() const {
603  checkInt(*this, "getInt()");
604  long int res = d_n->getNum();
605  FatalAssert(res >= INT_MIN && res <= INT_MAX,
606  "Rational::getInt(): arithmetic overflow on "+toString() +
607  OVERFLOW_MSG);
608  return (int)res;
609  }
610 
611  unsigned int Rational::getUnsigned() const {
612  checkInt(*this, "getUnsigned()");
613  long int res = d_n->getNum();
614  FatalAssert(res >= 0 && res <= (long int)UINT_MAX,
615  "Rational::getUnsigned(): arithmetic overflow on " + toString() +
616  OVERFLOW_MSG);
617  return (unsigned int)res;
618  }
619 
620 #ifdef _DEBUG_RATIONAL_
621  void Rational::printStats() {
622  int &num_created = getCreated();
623  int &num_deleted = getDeleted();
624  if(num_created % 1000 == 0 || num_deleted % 1000 == 0) {
625  std::cerr << "Rational(" << *d_n << "): created " << num_created
626  << ", deleted " << num_deleted
627  << ", currently alive " << num_created-num_deleted
628  << std::endl;
629  }
630  }
631 #endif
632 
633  bool operator==(const Rational &n1, const Rational &n2) {
634  return(*n1.d_n == *n2.d_n);
635  }
636 
637  bool operator<(const Rational &n1, const Rational &n2) {
638  return(*n1.d_n < *n2.d_n);
639  }
640 
641  bool operator<=(const Rational &n1, const Rational &n2) {
642  return(*n1.d_n <= *n2.d_n);
643  }
644 
645  bool operator>(const Rational &n1, const Rational &n2) {
646  return(*n1.d_n > *n2.d_n);
647  }
648 
649  bool operator>=(const Rational &n1, const Rational &n2) {
650  return(*n1.d_n >= *n2.d_n);
651  }
652 
653  bool operator!=(const Rational &n1, const Rational &n2) {
654  return(*n1.d_n != *n2.d_n);
655  }
656 
657  Rational operator+(const Rational &n1, const Rational &n2) {
658  return Rational(Rational::Impl(*n1.d_n + *n2.d_n));
659  }
660 
661  Rational operator-(const Rational &n1, const Rational &n2) {
662  return Rational(Rational::Impl((*n1.d_n) - (*n2.d_n)));
663  }
664 
665  Rational operator*(const Rational &n1, const Rational &n2) {
666  return Rational(Rational::Impl((*n1.d_n) * (*n2.d_n)));
667  }
668 
669  Rational operator/(const Rational &n1, const Rational &n2) {
670  return Rational(Rational::Impl(*n1.d_n / *n2.d_n));
671  }
672 
673  Rational operator%(const Rational &n1, const Rational &n2) {
674  return Rational(Rational::Impl(*n1.d_n % *n2.d_n));
675  }
676 
677  // Implementation of the forward-declared internal class
678  class Unsigned::Impl {
679  unsigned long d_num;
680  public:
681  //! Default Constructor
682  Impl(): d_num(0) { }
683  //! Copy constructor
684  Impl(const Impl &x) : d_num(x.d_num) { }
685  //! Constructor from an unsigned integer
686  Impl(unsigned long n): d_num(n) { }
687  //! Constructor from an int
688  Impl(int n): d_num(n) {
689  FatalAssert(n >= 0, "Attempt to create Unsigned from negative integer");
690  }
691 
692  //! Constructor from a string
693  Impl(const string &n, int base);
694  // Destructor
695  virtual ~Impl() { }
696  //! Get unsigned
697  unsigned long getUnsigned() const { return d_num; }
698 
699 
700  //! Equality
701  friend bool operator==(const Impl& x, const Impl& y) {
702  return (x.d_num == y.d_num);
703  }
704  //! Dis-equality
705  friend bool operator!=(const Impl& x, const Impl& y) {
706  return (x.d_num != y.d_num);
707  }
708  friend bool operator<(const Impl& x, const Impl& y) {
709  return x.d_num < y.d_num;
710  }
711 
712  friend bool operator<=(const Impl& x, const Impl& y) {
713  return (x.d_num <= y.d_num);
714  }
715 
716  friend bool operator>(const Impl& x, const Impl& y) {
717  return x.d_num > y.d_num;
718  }
719 
720  friend bool operator>=(const Impl& x, const Impl& y) {
721  return x.d_num >= y.d_num;
722  }
723 
724  friend Impl operator+(const Impl& x, const Impl& y) {
725  return Impl(uplus(x.d_num, y.d_num));
726  }
727 
728  friend Impl operator-(const Impl& x, const Impl& y) {
729  unsigned long n = unsigned_minus(x.d_num, y.d_num);
730  Impl res(n);
731  return res;
732  }
733 
734  friend Impl operator*(const Impl& x, const Impl& y) {
735  unsigned long n(umult(x.d_num, y.d_num));
736  return Impl(n);
737  }
738 
739  friend Impl operator/(const Impl& x, const Impl& y) {
740  DebugAssert(y.d_num != 0, "Impl::operator/: divisor is 0");
741  unsigned long n(x.d_num / y.d_num);
742  return Impl(n);
743  }
744 
745  friend Impl operator%(const Impl& x, const Impl& y) {
746  DebugAssert(y.d_num != 0,
747  "Impl % Impl: y must be non-zero");
748  return Impl(x.d_num % y.d_num);
749  }
750 
751  friend Impl operator<<(const Impl& x, unsigned y) {
752  unsigned long n(ushift(x.d_num, y));
753  return Impl(n);
754  }
755 
756  friend Impl operator&(const Impl& x, const Impl& y) {
757  return Impl(x.d_num & y.d_num);
758  }
759 
760  //! Print to string
761  string toString(int base = 10) const {
762  ostringstream ss;
763  if (d_num == 0) ss << "0";
764  else if (base == 10) {
765  ss << d_num;
766  }
767  else {
768  vector<int> vec;
769  long num = d_num;
770  while (num) {
771  vec.push_back(num % base);
772  num = num / base;
773  }
774  while (!vec.empty()) {
775  if (base > 10 && vec.back() > 10) {
776  ss << (char)('A' + (vec.back()-10));
777  }
778  else ss << vec.back();
779  vec.pop_back();
780  }
781  }
782  return(ss.str());
783  }
784 
785  //! Printing to ostream
786  friend ostream& operator<<(ostream& os, const Unsigned::Impl& n) {
787  return os << n.toString();
788  }
789 
790  };
791 
792  // Constructor from a pair of strings
793  Unsigned::Impl::Impl(const string &n, int base) {
794  d_num = strtoul(n.c_str(), NULL, base);
795  FatalAssert(d_num != ULONG_MAX,
796  "Unsigned::Impl(string): arithmetic overflow:"
797  "n = "+n+", base="+int2string(base)
798  +OVERFLOW_MSG);
799  }
800 
801  //! Default constructor
802  Unsigned::Unsigned() : d_n(new Impl) { }
803  //! Copy constructor
804  Unsigned::Unsigned(const Unsigned &n) : d_n(new Impl(*n.d_n)) { }
805 
806  //! Private constructor
807  Unsigned::Unsigned(const Impl& t): d_n(new Impl(t)) { }
808 
809  Unsigned::Unsigned(unsigned n): d_n(new Impl((unsigned long)n)) { }
810  Unsigned::Unsigned(int n): d_n(new Impl(n)) { }
811 
812  // Constructors from strings
813  Unsigned::Unsigned(const char* n, int base)
814  : d_n(new Impl(string(n), base)) { }
815  Unsigned::Unsigned(const string& n, int base)
816  : d_n(new Impl(n, base)) { }
817 
818  // Destructor
819  Unsigned::~Unsigned() {
820  delete d_n;
821  }
822 
823  // Assignment
824  Unsigned& Unsigned::operator=(const Unsigned& n) {
825  if(this == &n) return *this;
826  delete d_n;
827  d_n = new Impl(*n.d_n);
828  return *this;
829  }
830 
831  ostream &operator<<(ostream &os, const Unsigned &n) {
832  return(os << n.toString());
833  }
834 
835  Unsigned gcd(const Unsigned &x, const Unsigned &y) {
836  return Unsigned(Unsigned::Impl(ugcd(x.d_n->getUnsigned(), y.d_n->getUnsigned())));
837  }
838 
839  Unsigned gcd(const vector<Unsigned> &v) {
840  unsigned long g(1);
841  if(v.size() > 0) {
842  g = v[0].d_n->getUnsigned();
843  }
844  for(size_t i=1; i<v.size(); i++) {
845  if(g == 0) g = v[i].d_n->getUnsigned();
846  else if(v[i].d_n->getUnsigned() != 0)
847  g = ugcd(g, v[i].d_n->getUnsigned());
848  }
849  return Unsigned(Unsigned::Impl(g));
850  }
851 
852  Unsigned lcm(const Unsigned &x, const Unsigned &y) {
853  unsigned long g;
854  g = ulcm(x.d_n->getUnsigned(), y.d_n->getUnsigned());
855  return Unsigned(Unsigned::Impl(g));
856  }
857 
858  Unsigned lcm(const vector<Unsigned> &v) {
859  unsigned long g(1);
860  for(size_t i=0; i<v.size(); i++) {
861  if(v[i].d_n->getUnsigned() != 0)
862  g = ulcm(g, v[i].d_n->getUnsigned());
863  }
864  return Unsigned(Unsigned::Impl(g));
865  }
866 
867  Unsigned mod(const Unsigned &x, const Unsigned &y) {
868  unsigned long r = x.d_n->getUnsigned() % y.d_n->getUnsigned();
869  return(Unsigned(Unsigned::Impl(r)));
870  }
871 
872  Unsigned intRoot(const Unsigned& base, unsigned long int n) {
873  double b = base.d_n->getUnsigned();
874  double root = n;
875  root = 1/root;
876  b = ::pow(b, root);
877  unsigned long res = (unsigned long) ::floor(b);
878  if (::pow((long double)res, (int)n) == base.d_n->getUnsigned()) {
879  return Unsigned(Unsigned::Impl(res));
880  }
881  return Unsigned(Unsigned::Impl((unsigned long)0));
882  }
883 
884  string Unsigned::toString(int base) const {
885  return(d_n->toString(base));
886  }
887 
888  size_t Unsigned::hash() const {
890  return h(toString().c_str());
891  }
892 
893  void Unsigned::print() const {
894  cout << (*this) << endl;
895  }
896 
897  Unsigned &Unsigned::operator+=(const Unsigned &n2) {
898  *this = (*this) + n2;
899  return *this;
900  }
901 
902  Unsigned &Unsigned::operator-=(const Unsigned &n2) {
903  *this = (*this) - n2;
904  return *this;
905  }
906 
907  Unsigned &Unsigned::operator*=(const Unsigned &n2) {
908  *this = (*this) * n2;
909  return *this;
910  }
911 
912  Unsigned &Unsigned::operator/=(const Unsigned &n2) {
913  *this = (*this) / n2;
914  return *this;
915  }
916 
917  unsigned long Unsigned::getUnsigned() const {
918  return d_n->getUnsigned();
919  }
920 
921  bool operator==(const Unsigned &n1, const Unsigned &n2) {
922  return(*n1.d_n == *n2.d_n);
923  }
924 
925  bool operator<(const Unsigned &n1, const Unsigned &n2) {
926  return(*n1.d_n < *n2.d_n);
927  }
928 
929  bool operator<=(const Unsigned &n1, const Unsigned &n2) {
930  return(*n1.d_n <= *n2.d_n);
931  }
932 
933  bool operator>(const Unsigned &n1, const Unsigned &n2) {
934  return(*n1.d_n > *n2.d_n);
935  }
936 
937  bool operator>=(const Unsigned &n1, const Unsigned &n2) {
938  return(*n1.d_n >= *n2.d_n);
939  }
940 
941  bool operator!=(const Unsigned &n1, const Unsigned &n2) {
942  return(*n1.d_n != *n2.d_n);
943  }
944 
945  Unsigned operator+(const Unsigned &n1, const Unsigned &n2) {
946  return Unsigned(Unsigned::Impl(*n1.d_n + *n2.d_n));
947  }
948 
949  Unsigned operator-(const Unsigned &n1, const Unsigned &n2) {
950  return Unsigned(Unsigned::Impl((*n1.d_n) - (*n2.d_n)));
951  }
952 
953  Unsigned operator*(const Unsigned &n1, const Unsigned &n2) {
954  return Unsigned(Unsigned::Impl((*n1.d_n) * (*n2.d_n)));
955  }
956 
957  Unsigned operator/(const Unsigned &n1, const Unsigned &n2) {
958  return Unsigned(Unsigned::Impl(*n1.d_n / *n2.d_n));
959  }
960 
961  Unsigned operator%(const Unsigned &n1, const Unsigned &n2) {
962  return Unsigned(Unsigned::Impl(*n1.d_n % *n2.d_n));
963  }
964 
965  Unsigned operator<<(const Unsigned& n1, unsigned n2) {
966  return Unsigned(Unsigned::Impl(*n1.d_n << n2));
967  }
968 
969  Unsigned operator&(const Unsigned& n1, const Unsigned& n2) {
970  return Unsigned(Unsigned::Impl(*n1.d_n & *n2.d_n));
971  }
972 
973  Unsigned Rational::getUnsignedMP() const {
974  checkInt(*this, "getUnsignedMP()");
975  long int res = d_n->getNum();
976  FatalAssert(res >= 0 && res <= (long int)UINT_MAX,
977  "Rational::getUnsignedMP(): arithmetic overflow on " + toString() +
978  OVERFLOW_MSG);
979  return Unsigned((unsigned int)res);
980  }
981 
982 
983 } /* close namespace */
984 
985 #endif