00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include "globals.types.hh"
00024 #include "meta_programming.hh"
00025 #include "C_Integer.hh"
00026 #include <cassert>
00027
00028 namespace Parma_Polyhedra_Library {
00029
00030 namespace Checked {
00031
00032 template <typename T1, typename T2>
00033 struct Safe_Conversion : public False {
00034 };
00035 template <typename T>
00036 struct Safe_Conversion<T, T> : public True {
00037 };
00038
00039 #define safe_conversion(To, From) \
00040 template <> struct Safe_Conversion<To, From> : public True { }
00041
00042 safe_conversion(signed short, signed char);
00043 #if PPL_SIZEOF_CHAR < PPL_SIZEOF_SHORT
00044 safe_conversion(signed short, unsigned char);
00045 #endif
00046
00047 safe_conversion(signed int, signed char);
00048 safe_conversion(signed int, signed short);
00049 #if PPL_SIZEOF_CHAR < PPL_SIZEOF_INT
00050 safe_conversion(signed int, unsigned char);
00051 #endif
00052 #if PPL_SIZEOF_SHORT < PPL_SIZEOF_INT
00053 safe_conversion(signed int, unsigned short);
00054 #endif
00055
00056 safe_conversion(signed long, signed char);
00057 safe_conversion(signed long, signed short);
00058 safe_conversion(signed long, signed int);
00059 #if PPL_SIZEOF_CHAR < PPL_SIZEOF_LONG
00060 safe_conversion(signed long, unsigned char);
00061 #endif
00062 #if PPL_SIZEOF_SHORT < PPL_SIZEOF_LONG
00063 safe_conversion(signed long, unsigned short);
00064 #endif
00065 #if PPL_SIZEOF_INT < PPL_SIZEOF_LONG
00066 safe_conversion(signed long, unsigned int);
00067 #endif
00068
00069 safe_conversion(signed long long, signed char);
00070 safe_conversion(signed long long, signed short);
00071 safe_conversion(signed long long, signed int);
00072 safe_conversion(signed long long, signed long);
00073 #if PPL_SIZEOF_CHAR < PPL_SIZEOF_LONG_LONG
00074 safe_conversion(signed long long, unsigned char);
00075 #endif
00076 #if PPL_SIZEOF_SHORT < PPL_SIZEOF_LONG_LONG
00077 safe_conversion(signed long long, unsigned short);
00078 #endif
00079 #if PPL_SIZEOF_INT < PPL_SIZEOF_LONG_LONG
00080 safe_conversion(signed long long, unsigned int);
00081 #endif
00082 #if PPL_SIZEOF_LONG < PPL_SIZEOF_LONG_LONG
00083 safe_conversion(signed long long, unsigned long);
00084 #endif
00085
00086 safe_conversion(unsigned short, unsigned char);
00087
00088 safe_conversion(unsigned int, unsigned char);
00089 safe_conversion(unsigned int, unsigned short);
00090
00091 safe_conversion(unsigned long, unsigned char);
00092 safe_conversion(unsigned long, unsigned short);
00093 safe_conversion(unsigned long, unsigned int);
00094
00095 safe_conversion(unsigned long long, unsigned char);
00096 safe_conversion(unsigned long long, unsigned short);
00097 safe_conversion(unsigned long long, unsigned int);
00098 safe_conversion(unsigned long long, unsigned long);
00099
00100
00101 #if PPL_SIZEOF_CHAR <= PPL_SIZEOF_FLOAT - 2
00102 safe_conversion(float, signed char);
00103 safe_conversion(float, unsigned char);
00104 #endif
00105 #if PPL_SIZEOF_SHORT <= PPL_SIZEOF_FLOAT - 2
00106 safe_conversion(float, signed short);
00107 safe_conversion(float, unsigned short);
00108 #endif
00109 #if PPL_SIZEOF_INT <= PPL_SIZEOF_FLOAT - 2
00110 safe_conversion(float, signed int);
00111 safe_conversion(float, unsigned int);
00112 #endif
00113 #if PPL_SIZEOF_LONG <= PPL_SIZEOF_FLOAT - 2
00114 safe_conversion(float, signed long);
00115 safe_conversion(float, unsigned long);
00116 #endif
00117 #if PPL_SIZEOF_LONG_LONG <= PPL_SIZEOF_FLOAT - 2
00118 safe_conversion(float, signed long long);
00119 safe_conversion(float, unsigned long long);
00120 #endif
00121
00122 #if PPL_SIZEOF_CHAR <= PPL_SIZEOF_DOUBLE - 4
00123 safe_conversion(double, signed char);
00124 safe_conversion(double, unsigned char);
00125 #endif
00126 #if PPL_SIZEOF_SHORT <= PPL_SIZEOF_DOUBLE - 4
00127 safe_conversion(double, signed short);
00128 safe_conversion(double, unsigned short);
00129 #endif
00130 #if PPL_SIZEOF_INT <= PPL_SIZEOF_DOUBLE - 4
00131 safe_conversion(double, signed int);
00132 safe_conversion(double, unsigned int);
00133 #endif
00134 #if PPL_SIZEOF_LONG <= PPL_SIZEOF_DOUBLE - 4
00135 safe_conversion(double, signed long);
00136 safe_conversion(double, unsigned long);
00137 #endif
00138 #if PPL_SIZEOF_LONG_LONG <= PPL_SIZEOF_DOUBLE - 4
00139 safe_conversion(double, signed long long);
00140 safe_conversion(double, unsigned long long);
00141 #endif
00142 safe_conversion(double, float);
00143
00144 #if PPL_SIZEOF_CHAR <= PPL_SIZEOF_LONG_DOUBLE - 4
00145 safe_conversion(long double, signed char);
00146 safe_conversion(long double, unsigned char);
00147 #endif
00148 #if PPL_SIZEOF_SHORT <= PPL_SIZEOF_LONG_DOUBLE - 4
00149 safe_conversion(long double, signed short);
00150 safe_conversion(long double, unsigned short);
00151 #endif
00152 #if PPL_SIZEOF_INT <= PPL_SIZEOF_LONG_DOUBLE - 4
00153 safe_conversion(long double, signed int);
00154 safe_conversion(long double, unsigned int);
00155 #endif
00156 #if PPL_SIZEOF_LONG <= PPL_SIZEOF_LONG_DOUBLE - 4
00157 safe_conversion(long double, signed long);
00158 safe_conversion(long double, unsigned long);
00159 #endif
00160 #if PPL_SIZEOF_LONG_LONG <= PPL_SIZEOF_LONG_DOUBLE - 4
00161 safe_conversion(long double, signed long long);
00162 safe_conversion(long double, unsigned long long);
00163 #endif
00164 safe_conversion(long double, float);
00165 safe_conversion(long double, double);
00166
00167 safe_conversion(mpz_class, signed char);
00168 safe_conversion(mpz_class, signed short);
00169 safe_conversion(mpz_class, signed int);
00170 safe_conversion(mpz_class, signed long);
00171
00172 safe_conversion(mpz_class, unsigned char);
00173 safe_conversion(mpz_class, unsigned short);
00174 safe_conversion(mpz_class, unsigned int);
00175 safe_conversion(mpz_class, unsigned long);
00176
00177
00178 safe_conversion(mpq_class, signed char);
00179 safe_conversion(mpq_class, signed short);
00180 safe_conversion(mpq_class, signed int);
00181 safe_conversion(mpq_class, signed long);
00182
00183 safe_conversion(mpq_class, unsigned char);
00184 safe_conversion(mpq_class, unsigned short);
00185 safe_conversion(mpq_class, unsigned int);
00186 safe_conversion(mpq_class, unsigned long);
00187
00188 safe_conversion(mpq_class, float);
00189 safe_conversion(mpq_class, double);
00190
00191
00192 template <typename Policy, typename Type>
00193 struct FUNCTION_CLASS(construct)<Policy, Policy, Type, Type> {
00194 static inline Result function(Type& to, const Type& from, Rounding_Dir) {
00195 new (&to) Type(from);
00196 return V_EQ;
00197 }
00198 };
00199
00200 template <typename To_Policy, typename From_Policy, typename To, typename From>
00201 struct FUNCTION_CLASS(construct) {
00202 static inline Result function(To& to, const From& from, Rounding_Dir dir) {
00203 new (&to) To();
00204 return assign<To_Policy, From_Policy>(to, from, dir);
00205 }
00206 };
00207
00208 template <typename To_Policy, typename To>
00209 struct FUNCTION_CLASS(construct_special) {
00210 static inline Result function(To& to, Result r, Rounding_Dir dir) {
00211 new (&to) To();
00212 return assign_special<To_Policy>(to, r, dir);
00213 }
00214 };
00215
00216 template <typename To_Policy, typename From_Policy, typename To, typename From>
00217 inline Result
00218 assign_exact(To& to, const From& from, Rounding_Dir) {
00219 to = from;
00220 return V_EQ;
00221 }
00222
00223 template <typename To_Policy, typename From_Policy, typename Type>
00224 inline typename Enable_If<Is_Same<To_Policy, From_Policy>::value, void>::type
00225 copy_generic(Type& to, const Type& from) {
00226 to = from;
00227 }
00228
00229 template <typename To_Policy, typename From_Policy, typename To, typename From>
00230 inline Result
00231 abs_generic(To& to, const From& from, Rounding_Dir dir) {
00232 if (from < 0)
00233 return neg<To_Policy, From_Policy>(to, from, dir);
00234 else
00235 return assign<To_Policy, From_Policy>(to, from, dir);
00236 }
00237
00238 inline Result
00239 neg(Result r) {
00240 assert(!is_special(r));
00241 Result ret = static_cast<Result>(r & V_EQ);
00242 if (r & V_LT)
00243 ret = static_cast<Result>(ret | V_GT);
00244 if (r & V_GT)
00245 ret = static_cast<Result>(ret | V_LT);
00246 return ret;
00247 }
00248
00249 inline Result
00250 add(Result r1, Result r2) {
00251 assert(!is_special(r1));
00252 assert(!is_special(r2));
00253 if (r1 == V_EQ)
00254 return r2;
00255 if (r2 == V_EQ)
00256 return r1;
00257 if (((r1 & V_LT) && (r2 & V_GT))
00258 || ((r1 & V_GT) && (r2 & V_LT)))
00259 return V_LGE;
00260 return static_cast<Result>((((r1 & r2) & V_EQ) ? V_EQ : 0) |
00261 (r1 & (V_LT | V_GT)));
00262 }
00263
00264 inline Result
00265 sub(Result r1, Result r2) {
00266 return add(r1, neg(r2));
00267 }
00268
00269 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00270 typename To, typename From>
00271 inline void
00272 gcd_exact_noabs(To& to, const From& x, const From& y) {
00273 To nx = x;
00274 To ny = y;
00275 To rm;
00276 while (ny != 0) {
00277
00278
00279
00280 rem<To_Policy, From1_Policy, From2_Policy>(rm, nx, ny, ROUND_NOT_NEEDED);
00281 nx = ny;
00282 ny = rm;
00283 }
00284 to = nx;
00285 }
00286
00287 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00288 typename To, typename From1, typename From2>
00289 inline Result
00290 gcd_exact(To& to, const From1& x, const From2& y, Rounding_Dir dir) {
00291 gcd_exact_noabs<To_Policy, From1_Policy, From2_Policy>(to, x, y);
00292 return abs<To_Policy, To_Policy>(to, to, dir);
00293 }
00294
00295 template <typename To1_Policy, typename To2_Policy, typename To3_Policy,
00296 typename From1_Policy, typename From2_Policy,
00297 typename To1, typename To2, typename To3,
00298 typename From1, typename From2>
00299 inline Result
00300 gcdext_exact(To1& to, To2& s, To3& t, const From1& x, const From2& y,
00301 Rounding_Dir dir) {
00302
00303
00304
00305
00306 if (y == 0) {
00307 if (x == 0) {
00308 s = 0;
00309 t = 1;
00310 return V_EQ;
00311 }
00312 else {
00313 if (x < 0)
00314 s = -1;
00315 else
00316 s = 1;
00317 t = 0;
00318 return abs<To1_Policy, From1_Policy>(to, x, dir);
00319 }
00320 }
00321
00322 s = 1;
00323 t = 0;
00324 bool negative_x = x < 0;
00325 bool negative_y = y < 0;
00326
00327 Result r;
00328 r = abs<To1_Policy, From1_Policy>(to, x, dir);
00329 if (r != V_EQ)
00330 return r;
00331
00332 From2 ay;
00333 r = abs<To1_Policy, From2_Policy>(ay, y, dir);
00334 if (r != V_EQ)
00335 return r;
00336
00337
00338
00339
00340
00341 #define COPY_GMP
00342 #ifdef COPY_GMP
00343 if (to == ay)
00344 goto sign_check;
00345 #endif
00346
00347 {
00348 To2 v1 = 0;
00349 To3 v2 = 1;
00350 To1 v3 = static_cast<To1>(ay);
00351 while (true) {
00352 To1 q = to / v3;
00353
00354 To1 t3 = to - q*v3;
00355 To2 t1 = s - static_cast<To2>(q)*v1;
00356 To3 t2 = t - static_cast<To3>(q)*v2;
00357 s = v1;
00358 t = v2;
00359 to = v3;
00360 if (t3 == 0)
00361 break;
00362 v1 = t1;
00363 v2 = t2;
00364 v3 = t3;
00365 }
00366 }
00367
00368 #ifdef COPY_GMP
00369 sign_check:
00370 #endif
00371 if (negative_x) {
00372 r = neg<To2_Policy, To2_Policy>(s, s, dir);
00373 if (r != V_EQ)
00374 return r;
00375 }
00376 if (negative_y)
00377 return neg<To3_Policy, To3_Policy>(t, t, dir);
00378 return V_EQ;
00379 }
00380
00381 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00382 typename To, typename From1, typename From2>
00383 inline Result
00384 lcm_gcd_exact(To& to, const From1& x, const From2& y, Rounding_Dir dir) {
00385 if (x == 0 || y == 0) {
00386 to = 0;
00387 return V_EQ;
00388 }
00389 To nx, ny;
00390 Result r;
00391 r = abs<From1_Policy, From1_Policy>(nx, x, dir);
00392 if (r != V_EQ)
00393 return r;
00394 r = abs<From2_Policy, From2_Policy>(ny, y, dir);
00395 if (r != V_EQ)
00396 return r;
00397 To gcd;
00398 gcd_exact_noabs<To_Policy, From1_Policy, From2_Policy>(gcd, nx, ny);
00399
00400
00401
00402 div<To_Policy, From1_Policy, To_Policy>(to, nx, gcd, ROUND_NOT_NEEDED);
00403 return mul<To_Policy, To_Policy, From2_Policy>(to, to, ny, dir);
00404 }
00405
00406 template <typename Policy, typename Type>
00407 inline Result
00408 sgn_generic(const Type& x) {
00409 if (x > 0)
00410 return V_GT;
00411 if (x == 0)
00412 return V_EQ;
00413 return V_LT;
00414 }
00415
00416 template <typename T1, typename T2, typename Enable = void>
00417 struct Safe_Int_Comparison : public False {
00418 };
00419
00420 template <typename T1, typename T2>
00421 struct Safe_Int_Comparison<T1, T2, typename Enable_If<(C_Integer<T1>::value && C_Integer<T2>::value)>::type>
00422 : public Bool<(C_Integer<T1>::is_signed
00423 ? (C_Integer<T2>::is_signed
00424 || sizeof(T2) < sizeof(T1)
00425 || sizeof(T2) < sizeof(int))
00426 : (!C_Integer<T2>::is_signed
00427 || sizeof(T1) < sizeof(T2)
00428 || sizeof(T1) < sizeof(int)))> {
00429 };
00430
00431
00432 template <typename T1, typename T2>
00433 inline typename Enable_If<(Safe_Int_Comparison<T1, T2>::value
00434 || Safe_Conversion<T1, T2>::value
00435 || Safe_Conversion<T2, T1>::value), bool>::type
00436 lt(const T1& x, const T2& y) {
00437 return x < y;
00438 }
00439 template <typename T1, typename T2>
00440 inline typename Enable_If<(Safe_Int_Comparison<T1, T2>::value
00441 || Safe_Conversion<T1, T2>::value
00442 || Safe_Conversion<T2, T1>::value), bool>::type
00443 le(const T1& x, const T2& y) {
00444 return x <= y;
00445 }
00446 template <typename T1, typename T2>
00447 inline typename Enable_If<(Safe_Int_Comparison<T1, T2>::value
00448 || Safe_Conversion<T1, T2>::value
00449 || Safe_Conversion<T2, T1>::value), bool>::type
00450 eq(const T1& x, const T2& y) {
00451 return x == y;
00452 }
00453
00454 template <typename S, typename U>
00455 inline typename Enable_If<(!Safe_Int_Comparison<S, U>::value
00456 && C_Integer<U>::value
00457 && C_Integer<S>::is_signed), bool>::type
00458 lt(const S& x, const U& y) {
00459 return x < 0 || x < y;
00460 }
00461
00462 template <typename U, typename S>
00463 inline typename Enable_If<(!Safe_Int_Comparison<S, U>::value
00464 && C_Integer<U>::value
00465 && C_Integer<S>::is_signed), bool>::type
00466 lt(const U& x, const S& y) {
00467 return y >= 0 && x < y;
00468 }
00469
00470 template <typename S, typename U>
00471 inline typename Enable_If<(!Safe_Int_Comparison<S, U>::value
00472 && C_Integer<U>::value
00473 && C_Integer<S>::is_signed), bool>::type
00474 le(const S& x, const U& y) {
00475 return x < 0 || x <= y;
00476 }
00477
00478 template <typename U, typename S>
00479 inline typename Enable_If<(!Safe_Int_Comparison<S, U>::value
00480 && C_Integer<U>::value
00481 && C_Integer<S>::is_signed), bool>::type
00482 le(const U& x, const S& y) {
00483 return y >= 0 && x <= y;
00484 }
00485
00486 template <typename S, typename U>
00487 inline typename Enable_If<(!Safe_Int_Comparison<S, U>::value
00488 && C_Integer<U>::value
00489 && C_Integer<S>::is_signed), bool>::type
00490 eq(const S& x, const U& y) {
00491 return x >= 0 && x == y;
00492 }
00493
00494 template <typename U, typename S>
00495 inline typename Enable_If<(!Safe_Int_Comparison<S, U>::value
00496 && C_Integer<U>::value
00497 && C_Integer<S>::is_signed), bool>::type
00498 eq(const U& x, const S& y) {
00499 return y >= 0 && x == y;
00500 }
00501
00502 template <typename T1, typename T2>
00503 inline typename Enable_If<(!Safe_Conversion<T1, T2>::value
00504 && !Safe_Conversion<T2, T1>::value
00505 && (!C_Integer<T1>::value || !C_Integer<T2>::value)), bool>::type
00506 eq(const T1& x, const T2& y) {
00507 DIRTY_TEMP(T1, tmp);
00508 Result r = assign_r(tmp, y, ROUND_CHECK);
00509 return r == V_EQ && x == tmp;
00510 }
00511
00512 template <typename T1, typename T2>
00513 inline typename Enable_If<(!Safe_Conversion<T1, T2>::value
00514 && !Safe_Conversion<T2, T1>::value
00515 && (!C_Integer<T1>::value || !C_Integer<T2>::value)), bool>::type
00516 lt(const T1& x, const T2& y) {
00517 DIRTY_TEMP(T1, tmp);
00518 Result r = assign_r(tmp, y, ROUND_UP);
00519 switch (r) {
00520 case V_POS_OVERFLOW:
00521 case VC_PLUS_INFINITY:
00522 return true;
00523 case V_EQ:
00524 case V_LT:
00525 case V_LE:
00526 return x < tmp;
00527 default:
00528 return false;
00529 }
00530 }
00531
00532 template <typename T1, typename T2>
00533 inline typename Enable_If<(!Safe_Conversion<T1, T2>::value
00534 && !Safe_Conversion<T2, T1>::value
00535 && (!C_Integer<T1>::value || !C_Integer<T2>::value)), bool>::type
00536 le(const T1& x, const T2& y) {
00537 DIRTY_TEMP(T1, tmp);
00538 Result r = assign_r(tmp, y, static_cast<Rounding_Dir>(ROUND_UP | ROUND_FPU_CHECK_INEXACT));
00539 switch (r) {
00540 case V_POS_OVERFLOW:
00541 case VC_PLUS_INFINITY:
00542 return true;
00543 case V_EQ:
00544 return x <= tmp;
00545 case V_LT:
00546 return x < tmp;
00547 default:
00548 return false;
00549 }
00550 }
00551
00552 template <typename Policy1, typename Policy2,
00553 typename Type1, typename Type2>
00554 inline bool
00555 lt_p(const Type1& x, const Type2& y) {
00556 return lt(x, y);
00557 }
00558
00559 template <typename Policy1, typename Policy2,
00560 typename Type1, typename Type2>
00561 inline bool
00562 le_p(const Type1& x, const Type2& y) {
00563 return le(x, y);
00564 }
00565
00566 template <typename Policy1, typename Policy2,
00567 typename Type1, typename Type2>
00568 inline bool
00569 eq_p(const Type1& x, const Type2& y) {
00570 return eq(x, y);
00571 }
00572
00573 template <typename Policy1, typename Policy2,
00574 typename Type1, typename Type2>
00575 inline Result
00576 cmp_generic(const Type1& x, const Type2& y) {
00577 if (lt(y, x))
00578 return V_GT;
00579 if (lt(x, y))
00580 return V_LT;
00581 return V_EQ;
00582 }
00583
00584 template <typename Policy, typename Type>
00585 inline Result
00586 input_generic(Type& to, std::istream& is, Rounding_Dir dir) {
00587 DIRTY_TEMP0(mpq_class, q);
00588 Result r = input_mpq(q, is);
00589 if (is_special(r))
00590 return assign_special<Policy>(to, r, dir);
00591 if (r == V_EQ)
00592 return assign<Policy, void>(to, q, dir);
00593 assert(0);
00594 return VC_NAN;
00595 }
00596
00597 }
00598
00599 }