23 #ifdef RATIONAL_NATIVE
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."
41 static long int plus(
long int x,
long int y) {
43 FatalAssert(((x > 0) != (y > 0)) || ((x > 0) == (res > 0)),
44 "plus(x,y): arithmetic overflow" OVERFLOW_MSG);
49 static unsigned long uplus(
unsigned long x,
unsigned long y) {
50 unsigned long res = x+y;
52 "uplus(x,y): arithmetic overflow" OVERFLOW_MSG);
57 static unsigned long unsigned_minus(
unsigned long x,
unsigned long y) {
58 unsigned long res = x-y;
60 "unsigned_minus(x,y): arithmetic overflow" OVERFLOW_MSG);
65 static unsigned long umult(
unsigned long x,
unsigned long y) {
66 unsigned long res = x*y;
68 "umult(x,y): arithmetic overflow" OVERFLOW_MSG);
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);
78 "ushift(x,y): arithmetic overflow" OVERFLOW_MSG);
83 static unsigned long ugcd(
unsigned long n1,
unsigned long n2) {
101 static long int uminus(
long int x) {
102 FatalAssert(x == 0 || x != -x,
"uminus(x): arithmetic overflow"
108 static long int mult(
long int x,
long int y) {
110 FatalAssert(x==0 || res/x == y,
"mult(x,y): arithmetic overflow"
116 static long int gcd(
long int n1,
long int n2) {
120 if(n1 < 0) n1 = uminus(n1);
121 if(n2 < 0) n2 = uminus(n2);
130 long int r = n1 % n2;
138 static long int lcm(
long int n1,
long int n2) {
139 long int g = gcd(n1,n2);
140 return mult(n1/g, n2);
143 static long int ulcm(
unsigned long n1,
unsigned long n2) {
144 long int g = ugcd(n1,n2);
145 return umult(n1/g, n2);
149 class Rational::Impl {
156 Impl(): d_num(0), d_denom(1) { }
158 Impl(
const Impl &x) : d_num(x.d_num), d_denom(x.d_denom) { }
160 Impl(
unsigned long n) :d_num(n), d_denom(1) {
161 FatalAssert(d_num >= 0,
"Rational::Impl(unsigned long) - arithmetic overflow");
164 Impl(
long int n,
long int d): d_num(n), d_denom(d) { canonicalize(); }
166 Impl(
const string &n,
int base);
168 Impl(
const string &n,
const string& d,
int base);
172 long int getNum()
const {
return d_num; }
174 long int getDen()
const {
return d_denom; }
180 friend bool operator==(
const Impl& x,
const Impl& y) {
181 return (x.d_num == y.d_num && x.d_denom == y.d_denom);
184 friend bool operator!=(
const Impl& x,
const Impl& y) {
185 return (x.d_num != y.d_num || x.d_denom != y.d_denom);
191 friend bool operator<(
const Impl& x,
const Impl& y) {
193 return diff.d_num < 0;
196 friend bool operator<=(
const Impl& x,
const Impl& y) {
197 return (x == y || x < y);
200 friend bool operator>(
const Impl& x,
const Impl& y) {
202 return diff.d_num > 0;
205 friend bool operator>=(
const Impl& x,
const Impl& y) {
206 return (x == y || x > y);
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));
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)));
225 TRACE(
"rational",
" => ", res,
"");
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));
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);
254 friend Impl operator%(
const Impl& x,
const Impl& y) {
256 "Impl % Impl: x and y must be integers");
257 return Impl(x.getNum() % y.getNum(), 1);
261 string toString(
int base = 10)
const {
263 if (d_num == 0) ss <<
"0";
264 else if (base == 10) {
267 ss <<
"/" << d_denom;
273 vec.push_back(num % base);
276 while (!vec.empty()) {
277 if (base > 10 && vec.back() > 10) {
278 ss << (char)(
'A' + (vec.back()-10));
280 else ss << vec.back();
285 if (d_denom == 0) ss <<
"0";
289 vec.push_back(num % base);
292 while (!vec.empty()) {
293 if (base > 10 && vec.back() > 10) {
294 ss << (char)(
'A' + (vec.back()-10));
296 else ss << vec.back();
306 friend ostream&
operator<<(ostream& os,
const Rational::Impl& n) {
307 return os << n.toString();
313 void Rational::Impl::canonicalize() {
315 "Rational::Impl::canonicalize: bad denominator: "
321 d_num = uminus(d_num);
322 d_denom = uminus(d_denom);
324 long int d = gcd(d_num, d_denom);
333 Rational::Impl::Impl(
const string &n,
int base) {
335 for(i=0,iend=n.size(); i<iend && n[i] !=
'/'; ++i);
338 *
this = Impl(n.substr(0,i-1), n.substr(i+1,iend-1), base);
340 *
this = Impl(n,
"1", base);
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:"
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:"
359 return Impl(uminus(d_num), d_denom);
365 #ifdef _DEBUG_RATIONAL_
366 int &num_created = getCreated();
372 Rational::Rational(
const Rational &n) : d_n(new Impl(*n.d_n)) {
373 #ifdef _DEBUG_RATIONAL_
374 int &num_created = getCreated();
380 Rational::Rational(
const Unsigned &n) : d_n(new Impl(n.getUnsigned())) {
381 #ifdef _DEBUG_RATIONAL_
382 int &num_created = getCreated();
389 Rational::Rational(
const Impl& t): d_n(new Impl(t)) {
390 #ifdef _DEBUG_RATIONAL_
391 int &num_created = getCreated();
396 Rational::Rational(
int n,
int d): d_n(new Impl(n, d)) {
397 #ifdef _DEBUG_RATIONAL_
398 int &num_created = getCreated();
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();
412 Rational::Rational(
const string& n,
int base)
413 : d_n(new Impl(n, base)) {
414 #ifdef _DEBUG_RATIONAL_
415 int &num_created = getCreated();
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();
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();
437 Rational::~Rational() {
439 #ifdef _DEBUG_RATIONAL_
440 int &num_deleted = getDeleted();
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)); }
452 bool Rational::isInteger()
const {
return 1 == d_n->getDen(); }
455 Rational& Rational::operator=(
const Rational& n) {
456 if(
this == &n)
return *
this;
458 d_n =
new Impl(*n.d_n);
462 ostream &
operator<<(ostream &os,
const Rational &n) {
463 return(os << n.toString());
469 static void checkInt(
const Rational& n,
const string& funName) {
471 (
"CVC3::Rational::" + funName
472 +
": argument is not an integer: " + n.toString()).c_str());
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));
485 Rational gcd(
const vector<Rational> &v) {
488 checkInt(v[0],
"gcd(vector<Rational>[0])");
489 g = v[0].d_n->getNum();
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());
497 return Rational(Rational::Impl(g,1));
500 Rational lcm(
const Rational &x,
const Rational &y) {
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));
508 Rational lcm(
const vector<Rational> &v) {
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());
515 return Rational(Rational::Impl(g,1));
518 Rational
abs(
const Rational &x) {
519 long int n(x.d_n->getNum());
521 return Rational(Rational::Impl(-n, x.d_n->getDen()));
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));
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);
540 return Rational(Rational::Impl(q,1));
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)));
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();
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));
561 return Rational(Rational::Impl((
long)0,(
long)1));
564 string Rational::toString(
int base)
const {
565 return(d_n->toString(base));
568 size_t Rational::hash()
const {
570 return h(toString().c_str());
573 void Rational::print()
const {
574 cout << (*this) <<
endl;
579 return Rational(Rational::Impl(-(d_n->getNum()), d_n->getDen()));
582 Rational &Rational::operator+=(
const Rational &n2) {
583 *
this = (*this) + n2;
587 Rational &Rational::operator-=(
const Rational &n2) {
588 *
this = (*this) - n2;
592 Rational &Rational::operator*=(
const Rational &n2) {
593 *
this = (*this) * n2;
597 Rational &Rational::operator/=(
const Rational &n2) {
598 *
this = (*this) / n2;
602 int Rational::getInt()
const {
603 checkInt(*
this,
"getInt()");
604 long int res = d_n->getNum();
606 "Rational::getInt(): arithmetic overflow on "+toString() +
611 unsigned int Rational::getUnsigned()
const {
612 checkInt(*
this,
"getUnsigned()");
613 long int res = d_n->getNum();
615 "Rational::getUnsigned(): arithmetic overflow on " + toString() +
617 return (
unsigned int)res;
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
633 bool operator==(
const Rational &n1,
const Rational &n2) {
634 return(*n1.d_n == *n2.d_n);
637 bool operator<(
const Rational &n1,
const Rational &n2) {
638 return(*n1.d_n < *n2.d_n);
641 bool operator<=(
const Rational &n1,
const Rational &n2) {
642 return(*n1.d_n <= *n2.d_n);
645 bool operator>(
const Rational &n1,
const Rational &n2) {
646 return(*n1.d_n > *n2.d_n);
649 bool operator>=(
const Rational &n1,
const Rational &n2) {
650 return(*n1.d_n >= *n2.d_n);
653 bool operator!=(
const Rational &n1,
const Rational &n2) {
654 return(*n1.d_n != *n2.d_n);
657 Rational
operator+(
const Rational &n1,
const Rational &n2) {
658 return Rational(Rational::Impl(*n1.d_n + *n2.d_n));
661 Rational
operator-(
const Rational &n1,
const Rational &n2) {
662 return Rational(Rational::Impl((*n1.d_n) - (*n2.d_n)));
665 Rational
operator*(
const Rational &n1,
const Rational &n2) {
666 return Rational(Rational::Impl((*n1.d_n) * (*n2.d_n)));
669 Rational
operator/(
const Rational &n1,
const Rational &n2) {
670 return Rational(Rational::Impl(*n1.d_n / *n2.d_n));
673 Rational operator%(
const Rational &n1,
const Rational &n2) {
674 return Rational(Rational::Impl(*n1.d_n % *n2.d_n));
678 class Unsigned::Impl {
684 Impl(
const Impl &x) : d_num(x.d_num) { }
686 Impl(
unsigned long n): d_num(n) { }
688 Impl(
int n): d_num(n) {
689 FatalAssert(n >= 0,
"Attempt to create Unsigned from negative integer");
693 Impl(
const string &n,
int base);
697 unsigned long getUnsigned()
const {
return d_num; }
701 friend bool operator==(
const Impl& x,
const Impl& y) {
702 return (x.d_num == y.d_num);
705 friend bool operator!=(
const Impl& x,
const Impl& y) {
706 return (x.d_num != y.d_num);
708 friend bool operator<(
const Impl& x,
const Impl& y) {
709 return x.d_num < y.d_num;
712 friend bool operator<=(
const Impl& x,
const Impl& y) {
713 return (x.d_num <= y.d_num);
716 friend bool operator>(
const Impl& x,
const Impl& y) {
717 return x.d_num > y.d_num;
720 friend bool operator>=(
const Impl& x,
const Impl& y) {
721 return x.d_num >= y.d_num;
724 friend Impl
operator+(
const Impl& x,
const Impl& y) {
725 return Impl(uplus(x.d_num, y.d_num));
728 friend Impl
operator-(
const Impl& x,
const Impl& y) {
729 unsigned long n = unsigned_minus(x.d_num, y.d_num);
734 friend Impl
operator*(
const Impl& x,
const Impl& y) {
735 unsigned long n(umult(x.d_num, y.d_num));
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);
745 friend Impl operator%(
const Impl& x,
const Impl& y) {
747 "Impl % Impl: y must be non-zero");
748 return Impl(x.d_num % y.d_num);
751 friend Impl
operator<<(
const Impl& x,
unsigned y) {
752 unsigned long n(ushift(x.d_num, y));
756 friend Impl operator&(
const Impl& x,
const Impl& y) {
757 return Impl(x.d_num & y.d_num);
761 string toString(
int base = 10)
const {
763 if (d_num == 0) ss <<
"0";
764 else if (base == 10) {
771 vec.push_back(num % base);
774 while (!vec.empty()) {
775 if (base > 10 && vec.back() > 10) {
776 ss << (char)(
'A' + (vec.back()-10));
778 else ss << vec.back();
786 friend ostream&
operator<<(ostream& os,
const Unsigned::Impl& n) {
787 return os << n.toString();
793 Unsigned::Impl::Impl(
const string &n,
int base) {
794 d_num = strtoul(n.c_str(), NULL, base);
796 "Unsigned::Impl(string): arithmetic overflow:"
802 Unsigned::Unsigned() : d_n(new Impl) { }
804 Unsigned::Unsigned(
const Unsigned &n) : d_n(new Impl(*n.d_n)) { }
807 Unsigned::Unsigned(
const Impl& t): d_n(new Impl(t)) { }
809 Unsigned::Unsigned(
unsigned n): d_n(new Impl((unsigned long)n)) { }
810 Unsigned::Unsigned(
int n): d_n(new Impl(n)) { }
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)) { }
819 Unsigned::~Unsigned() {
824 Unsigned& Unsigned::operator=(
const Unsigned& n) {
825 if(
this == &n)
return *
this;
827 d_n =
new Impl(*n.d_n);
831 ostream &
operator<<(ostream &os,
const Unsigned &n) {
832 return(os << n.toString());
835 Unsigned gcd(
const Unsigned &x,
const Unsigned &y) {
836 return Unsigned(Unsigned::Impl(ugcd(x.d_n->getUnsigned(), y.d_n->getUnsigned())));
839 Unsigned gcd(
const vector<Unsigned> &v) {
842 g = v[0].d_n->getUnsigned();
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());
849 return Unsigned(Unsigned::Impl(g));
852 Unsigned lcm(
const Unsigned &x,
const Unsigned &y) {
854 g = ulcm(x.d_n->getUnsigned(), y.d_n->getUnsigned());
855 return Unsigned(Unsigned::Impl(g));
858 Unsigned lcm(
const vector<Unsigned> &v) {
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());
864 return Unsigned(Unsigned::Impl(g));
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)));
872 Unsigned intRoot(
const Unsigned& base,
unsigned long int n) {
873 double b = base.d_n->getUnsigned();
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));
881 return Unsigned(Unsigned::Impl((
unsigned long)0));
884 string Unsigned::toString(
int base)
const {
885 return(d_n->toString(base));
888 size_t Unsigned::hash()
const {
890 return h(toString().c_str());
893 void Unsigned::print()
const {
894 cout << (*this) <<
endl;
897 Unsigned &Unsigned::operator+=(
const Unsigned &n2) {
898 *
this = (*this) + n2;
902 Unsigned &Unsigned::operator-=(
const Unsigned &n2) {
903 *
this = (*this) - n2;
907 Unsigned &Unsigned::operator*=(
const Unsigned &n2) {
908 *
this = (*this) * n2;
912 Unsigned &Unsigned::operator/=(
const Unsigned &n2) {
913 *
this = (*this) / n2;
917 unsigned long Unsigned::getUnsigned()
const {
918 return d_n->getUnsigned();
921 bool operator==(
const Unsigned &n1,
const Unsigned &n2) {
922 return(*n1.d_n == *n2.d_n);
925 bool operator<(
const Unsigned &n1,
const Unsigned &n2) {
926 return(*n1.d_n < *n2.d_n);
929 bool operator<=(
const Unsigned &n1,
const Unsigned &n2) {
930 return(*n1.d_n <= *n2.d_n);
933 bool operator>(
const Unsigned &n1,
const Unsigned &n2) {
934 return(*n1.d_n > *n2.d_n);
937 bool operator>=(
const Unsigned &n1,
const Unsigned &n2) {
938 return(*n1.d_n >= *n2.d_n);
941 bool operator!=(
const Unsigned &n1,
const Unsigned &n2) {
942 return(*n1.d_n != *n2.d_n);
945 Unsigned
operator+(
const Unsigned &n1,
const Unsigned &n2) {
946 return Unsigned(Unsigned::Impl(*n1.d_n + *n2.d_n));
949 Unsigned
operator-(
const Unsigned &n1,
const Unsigned &n2) {
950 return Unsigned(Unsigned::Impl((*n1.d_n) - (*n2.d_n)));
953 Unsigned
operator*(
const Unsigned &n1,
const Unsigned &n2) {
954 return Unsigned(Unsigned::Impl((*n1.d_n) * (*n2.d_n)));
957 Unsigned
operator/(
const Unsigned &n1,
const Unsigned &n2) {
958 return Unsigned(Unsigned::Impl(*n1.d_n / *n2.d_n));
961 Unsigned operator%(
const Unsigned &n1,
const Unsigned &n2) {
962 return Unsigned(Unsigned::Impl(*n1.d_n % *n2.d_n));
965 Unsigned
operator<<(
const Unsigned& n1,
unsigned n2) {
966 return Unsigned(Unsigned::Impl(*n1.d_n << n2));
969 Unsigned operator&(
const Unsigned& n1,
const Unsigned& n2) {
970 return Unsigned(Unsigned::Impl(*n1.d_n & *n2.d_n));
973 Unsigned Rational::getUnsignedMP()
const {
974 checkInt(*
this,
"getUnsignedMP()");
975 long int res = d_n->getNum();
977 "Rational::getUnsignedMP(): arithmetic overflow on " + toString() +
979 return Unsigned((
unsigned int)res);