00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef PPL_checked_float_inlines_hh
00024 #define PPL_checked_float_inlines_hh 1
00025
00026 #ifndef __alpha
00027 #include <cmath>
00028 #endif
00029
00030 namespace Parma_Polyhedra_Library {
00031
00032 namespace Checked {
00033
00034 inline float
00035 fma(float x, float y, float z) {
00036 #if PPL_HAVE_DECL_FMAF && !defined(__alpha)
00037 return ::fmaf(x, y, z);
00038 #else
00039 return x*y + z;
00040 #endif
00041 }
00042
00043 inline double
00044 fma(double x, double y, double z) {
00045 #if PPL_HAVE_DECL_FMA && !defined(__alpha)
00046 return ::fma(x, y, z);
00047 #else
00048 return x*y + z;
00049 #endif
00050 }
00051
00052 inline long double
00053 fma(long double x, long double y, long double z) {
00054 #if PPL_HAVE_DECL_FMAL && !defined(__alpha)
00055 return ::fmal(x, y, z);
00056 #else
00057 return x*y + z;
00058 #endif
00059 }
00060
00061 #if PPL_HAVE_DECL_RINTF
00062 inline float
00063 rint(float x) {
00064 return ::rintf(x);
00065 }
00066 #endif
00067
00068 inline double
00069 rint(double x) {
00070 return ::rint(x);
00071 }
00072
00073 #if PPL_HAVE_DECL_RINTL
00074 inline long double
00075 rint(long double x) {
00076 return ::rintl(x);
00077 }
00078 #elif !PPL_CXX_PROVIDES_PROPER_LONG_DOUBLE
00079
00080
00081 inline long double
00082 rint(long double x) {
00083 return ::rint(x);
00084 }
00085 #elif defined(__i386__) && (defined(__GNUC__) || defined(__INTEL_COMPILER))
00086
00087
00088 inline long double
00089 rint(long double x) {
00090 long double i;
00091 __asm__ ("frndint" : "=t" (i) : "0" (x));
00092 return i;
00093 }
00094 #endif
00095
00096 inline bool
00097 fpu_direct_rounding(Rounding_Dir dir) {
00098 return round_direct(dir) || round_ignore(dir);
00099 }
00100
00101 inline bool
00102 fpu_inverse_rounding(Rounding_Dir dir) {
00103 return round_inverse(dir);
00104 }
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 inline float
00131 limit_precision(float v) {
00132 float x = v;
00133 avoid_cse(x);
00134 return x;
00135 }
00136
00137 inline double
00138 limit_precision(double v) {
00139 double x = v;
00140 avoid_cse(x);
00141 return x;
00142 }
00143
00144 inline long double
00145 limit_precision(long double v) {
00146 return v;
00147 }
00148
00149 template <typename Policy, typename T>
00150 inline Result
00151 classify_float(const T v, bool nan, bool inf, bool sign) {
00152 Float<T> f(v);
00153 if ((nan || sign) && CHECK_P(Policy::has_nan, f.u.binary.is_nan()))
00154 return VC_NAN;
00155 if (inf) {
00156 int i = CHECK_P(Policy::has_infinity, f.u.binary.is_inf());
00157 if (i < 0)
00158 return VC_MINUS_INFINITY;
00159 if (i > 0)
00160 return VC_PLUS_INFINITY;
00161 }
00162 if (sign) {
00163 if (v < 0)
00164 return V_LT;
00165 if (v > 0)
00166 return V_GT;
00167 return V_EQ;
00168 }
00169 return VC_NORMAL;
00170 }
00171
00172 template <typename Policy, typename T>
00173 inline bool
00174 is_nan_float(const T v) {
00175 Float<T> f(v);
00176 return CHECK_P(Policy::has_nan, f.u.binary.is_nan());
00177 }
00178
00179 template <typename Policy, typename T>
00180 inline int
00181 is_inf_float(const T v) {
00182 Float<T> f(v);
00183 return CHECK_P(Policy::has_infinity, f.u.binary.is_inf());
00184 }
00185 template <typename Policy, typename T>
00186 inline bool
00187 is_minf_float(const T v) {
00188 return is_inf_float<Policy>(v) < 0;
00189 }
00190
00191 template <typename Policy, typename T>
00192 inline bool
00193 is_pinf_float(const T v) {
00194 return is_inf_float<Policy>(v) > 0;
00195 }
00196
00197
00198 template <typename Policy, typename T>
00199 inline bool
00200 is_int_float(const T v) {
00201 return rint(v) == v;
00202 }
00203
00204 template <typename Policy, typename T>
00205 inline Result
00206 assign_special_float(T& v, Result r, Rounding_Dir) {
00207 switch (classify(r)) {
00208 case VC_MINUS_INFINITY:
00209 v = -HUGE_VAL;
00210 break;
00211 case VC_PLUS_INFINITY:
00212 v = HUGE_VAL;
00213 break;
00214 case VC_NAN:
00215 v = NAN;
00216 return r;
00217 default:
00218 break;
00219 }
00220 return V_EQ;
00221 }
00222
00223 template <typename T>
00224 inline void
00225 pred_float(T& v) {
00226 Float<T> f(v);
00227 assert(!f.u.binary.is_nan());
00228 assert(f.u.binary.is_inf() >= 0);
00229 if (f.u.binary.is_zero() > 0) {
00230 f.u.binary.negate();
00231 f.u.binary.inc();
00232 }
00233 else if (f.u.binary.sign_bit()) {
00234 f.u.binary.inc();
00235 }
00236 else {
00237 f.u.binary.dec();
00238 }
00239 v = f.value();
00240 }
00241
00242 template <typename T>
00243 inline void
00244 succ_float(T& v) {
00245 Float<T> f(v);
00246 assert(!f.u.binary.is_nan());
00247 assert(f.u.binary.is_inf() <= 0);
00248 if (f.u.binary.is_zero() < 0) {
00249 f.u.binary.negate();
00250 f.u.binary.inc();
00251 }
00252 else if (!f.u.binary.sign_bit()) {
00253 f.u.binary.inc();
00254 }
00255 else {
00256 f.u.binary.dec();
00257 }
00258 v = f.value();
00259 }
00260
00261 template <typename Policy, typename To>
00262 inline Result
00263 round_lt_float(To& to, Rounding_Dir dir) {
00264 if (round_down(dir)) {
00265 pred_float(to);
00266 return V_GT;
00267 }
00268 return V_LT;
00269 }
00270
00271 template <typename Policy, typename To>
00272 inline Result
00273 round_gt_float(To& to, Rounding_Dir dir) {
00274 if (round_up(dir)) {
00275 succ_float(to);
00276 return V_LT;
00277 }
00278 return V_GT;
00279 }
00280
00281
00282 template <typename Policy>
00283 inline void
00284 prepare_inexact(Rounding_Dir dir) {
00285 if (Policy::fpu_check_inexact && round_fpu_check_inexact(dir))
00286 fpu_reset_inexact();
00287 }
00288
00289 template <typename Policy>
00290 inline Result
00291 result_relation(Rounding_Dir dir) {
00292 if (Policy::fpu_check_inexact && round_fpu_check_inexact(dir)) {
00293 switch (fpu_check_inexact()) {
00294 case 0:
00295 return V_EQ;
00296 case -1:
00297 goto unknown;
00298 case 1:
00299 break;
00300 }
00301 switch (round_dir(dir)) {
00302 case ROUND_DOWN:
00303 return V_GT;
00304 case ROUND_UP:
00305 return V_LT;
00306 default:
00307 return V_NE;
00308 }
00309 }
00310 else {
00311 unknown:
00312 switch (round_dir(dir)) {
00313 case ROUND_DOWN:
00314 return V_GE;
00315 case ROUND_UP:
00316 return V_LE;
00317 default:
00318 return V_LGE;
00319 }
00320 }
00321 }
00322
00323 template <typename To_Policy, typename From_Policy, typename To, typename From>
00324 inline Result
00325 assign_float_float_exact(To& to, const From from, Rounding_Dir) {
00326 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00327 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00328 to = from;
00329 return V_EQ;
00330 }
00331
00332 template <typename To_Policy, typename From_Policy, typename To, typename From>
00333 inline Result
00334 assign_float_float_inexact(To& to, const From from, Rounding_Dir dir) {
00335 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00336 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00337 prepare_inexact<To_Policy>(dir);
00338 if (fpu_direct_rounding(dir))
00339 to = from;
00340 else if (fpu_inverse_rounding(dir))
00341 to = -static_cast<To>(-from);
00342 else {
00343 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00344 avoid_cse(from);
00345 to = from;
00346 avoid_cse(to);
00347 fpu_restore_rounding_direction(old);
00348 }
00349 return result_relation<To_Policy>(dir);
00350 }
00351
00352 template <typename To_Policy, typename From_Policy, typename To, typename From>
00353 inline Result
00354 assign_float_float(To& to, const From from, Rounding_Dir dir) {
00355 if (sizeof(From) > sizeof(To))
00356 return assign_float_float_inexact<To_Policy, From_Policy>(to, from, dir);
00357 else
00358 return assign_float_float_exact<To_Policy, From_Policy>(to, from, dir);
00359 }
00360
00361 template <typename To_Policy, typename From_Policy, typename Type>
00362 inline Result
00363 floor_float(Type& to, const Type from, Rounding_Dir) {
00364 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00365 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00366 if (fpu_direct_rounding(ROUND_DOWN))
00367 to = rint(from);
00368 else if (fpu_inverse_rounding(ROUND_DOWN))
00369 to = -limit_precision(rint(-from));
00370 else {
00371 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(ROUND_DOWN));
00372 avoid_cse(from);
00373 to = rint(from);
00374 avoid_cse(to);
00375 fpu_restore_rounding_direction(old);
00376 }
00377 return V_EQ;
00378 }
00379
00380 template <typename To_Policy, typename From_Policy, typename Type>
00381 inline Result
00382 ceil_float(Type& to, const Type from, Rounding_Dir) {
00383 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00384 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00385 if (fpu_direct_rounding(ROUND_UP))
00386 to = rint(from);
00387 else if (fpu_inverse_rounding(ROUND_UP))
00388 to = -limit_precision(rint(-from));
00389 else {
00390 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(ROUND_UP));
00391 avoid_cse(from);
00392 to = rint(from);
00393 avoid_cse(to);
00394 fpu_restore_rounding_direction(old);
00395 }
00396 return V_EQ;
00397 }
00398
00399 template <typename To_Policy, typename From_Policy, typename Type>
00400 inline Result
00401 trunc_float(Type& to, const Type from, Rounding_Dir dir) {
00402 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00403 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00404 if (from >= 0)
00405 return floor<To_Policy, From_Policy>(to, from, dir);
00406 else
00407 return ceil<To_Policy, From_Policy>(to, from, dir);
00408 }
00409
00410 template <typename To_Policy, typename From_Policy, typename Type>
00411 inline Result
00412 neg_float(Type& to, const Type from, Rounding_Dir) {
00413 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00414 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00415 to = -from;
00416 return V_EQ;
00417 }
00418
00419 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00420 typename Type>
00421 inline Result
00422 add_float(Type& to, const Type x, const Type y, Rounding_Dir dir) {
00423 if (To_Policy::check_inf_add_inf
00424 && is_inf_float<From1_Policy>(x) && x == -y)
00425 return assign_special<To_Policy>(to, V_INF_ADD_INF, ROUND_IGNORE);
00426 prepare_inexact<To_Policy>(dir);
00427 if (fpu_direct_rounding(dir))
00428 to = x + y;
00429 else if (fpu_inverse_rounding(dir))
00430 to = -limit_precision(-x - y);
00431 else {
00432 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00433 avoid_cse(x);
00434 avoid_cse(y);
00435 to = x + y;
00436 avoid_cse(to);
00437 fpu_restore_rounding_direction(old);
00438 }
00439 if (To_Policy::check_nan_result && is_nan<To_Policy>(to))
00440 return VC_NAN;
00441 return result_relation<To_Policy>(dir);
00442 }
00443
00444 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00445 typename Type>
00446 inline Result
00447 sub_float(Type& to, const Type x, const Type y, Rounding_Dir dir) {
00448 if (To_Policy::check_inf_sub_inf
00449 && is_inf_float<From1_Policy>(x) && x == y)
00450 return assign_special<To_Policy>(to, V_INF_SUB_INF, ROUND_IGNORE);
00451 prepare_inexact<To_Policy>(dir);
00452 if (fpu_direct_rounding(dir))
00453 to = x - y;
00454 else if (fpu_inverse_rounding(dir))
00455 to = -limit_precision(y - x);
00456 else {
00457 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00458 avoid_cse(x);
00459 avoid_cse(y);
00460 to = x - y;
00461 avoid_cse(to);
00462 fpu_restore_rounding_direction(old);
00463 }
00464 if (To_Policy::check_nan_result && is_nan<To_Policy>(to))
00465 return VC_NAN;
00466 return result_relation<To_Policy>(dir);
00467 }
00468
00469 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00470 typename Type>
00471 inline Result
00472 mul_float(Type& to, const Type x, const Type y, Rounding_Dir dir) {
00473 if (To_Policy::check_inf_mul_zero
00474 && ((x == 0 && is_inf_float<From2_Policy>(y)) ||
00475 (y == 0 && is_inf_float<From1_Policy>(x))))
00476 return assign_special<To_Policy>(to, V_INF_MUL_ZERO, ROUND_IGNORE);
00477 prepare_inexact<To_Policy>(dir);
00478 if (fpu_direct_rounding(dir))
00479 to = x * y;
00480 else if (fpu_inverse_rounding(dir))
00481 to = -limit_precision(x * -y);
00482 else {
00483 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00484 avoid_cse(x);
00485 avoid_cse(y);
00486 to = x * y;
00487 avoid_cse(to);
00488 fpu_restore_rounding_direction(old);
00489 }
00490 if (To_Policy::check_nan_result && is_nan<To_Policy>(to))
00491 return VC_NAN;
00492 return result_relation<To_Policy>(dir);
00493 }
00494
00495 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00496 typename Type>
00497 inline Result
00498 div_float(Type& to, const Type x, const Type y, Rounding_Dir dir) {
00499 if (To_Policy::check_inf_div_inf
00500 && is_inf_float<From1_Policy>(x) && is_inf_float<From2_Policy>(y))
00501 return assign_special<To_Policy>(to, V_INF_DIV_INF, ROUND_IGNORE);
00502 if (To_Policy::check_div_zero && y == 0)
00503 return assign_special<To_Policy>(to, V_DIV_ZERO, ROUND_IGNORE);
00504 prepare_inexact<To_Policy>(dir);
00505 if (fpu_direct_rounding(dir))
00506 to = x / y;
00507 else if (fpu_inverse_rounding(dir))
00508 to = -limit_precision(x / -y);
00509 else {
00510 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00511 avoid_cse(x);
00512 avoid_cse(y);
00513 to = x / y;
00514 avoid_cse(to);
00515 fpu_restore_rounding_direction(old);
00516 }
00517 if (To_Policy::check_nan_result && is_nan<To_Policy>(to))
00518 return VC_NAN;
00519 return result_relation<To_Policy>(dir);
00520 }
00521
00522 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00523 typename Type>
00524 inline Result
00525 idiv_float(Type& to, const Type x, const Type y, Rounding_Dir dir) {
00526 Type temp;
00527
00528 dir = round_dir(dir);
00529 Result r = div<To_Policy, From1_Policy, From2_Policy>(temp, x, y, dir);
00530 if (is_special(r)) {
00531 to = temp;
00532 return r;
00533 }
00534 Result r1 = trunc<To_Policy, To_Policy>(to, temp, ROUND_NOT_NEEDED);
00535 assert(r1 == V_EQ);
00536 if (r == V_EQ || to != temp)
00537 return r1;
00538
00539 return dir == ROUND_UP ? V_LE : V_GE;
00540 }
00541
00542 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00543 typename Type>
00544 inline Result
00545 rem_float(Type& to, const Type x, const Type y, Rounding_Dir) {
00546 if (To_Policy::check_inf_mod && is_inf_float<From1_Policy>(x))
00547 return assign_special<To_Policy>(to, V_INF_MOD, ROUND_IGNORE);
00548 if (To_Policy::check_div_zero && y == 0)
00549 return assign_special<To_Policy>(to, V_MOD_ZERO, ROUND_IGNORE);
00550 to = std::fmod(x, y);
00551 if (To_Policy::check_nan_result && is_nan<To_Policy>(to))
00552 return VC_NAN;
00553 return V_EQ;
00554 }
00555
00556 struct Float_2exp {
00557 const_bool_nodef(has_nan, false);
00558 const_bool_nodef(has_infinity, false);
00559 };
00560
00561 template <typename To_Policy, typename From_Policy, typename Type>
00562 inline Result
00563 mul2exp_float(Type& to, const Type x, int exp, Rounding_Dir dir) {
00564 if (To_Policy::check_nan_result && is_nan<From_Policy>(x))
00565 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00566 if (exp < 0)
00567 return div2exp<To_Policy, From_Policy>(to, x, -exp, dir);
00568 assert(static_cast<unsigned int>(exp) < sizeof(unsigned long long) * 8);
00569 return mul<To_Policy, From_Policy, Float_2exp>(to, x, static_cast<Type>(1ULL << exp), dir);
00570 }
00571
00572 template <typename To_Policy, typename From_Policy, typename Type>
00573 inline Result
00574 div2exp_float(Type& to, const Type x, int exp, Rounding_Dir dir) {
00575 if (To_Policy::check_nan_result && is_nan<From_Policy>(x))
00576 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00577 if (exp < 0)
00578 return mul2exp<To_Policy, From_Policy>(to, x, -exp, dir);
00579 assert(static_cast<unsigned int>(exp) < sizeof(unsigned long long) * 8);
00580 return div<To_Policy, From_Policy, Float_2exp>(to, x, static_cast<Type>(1ULL << exp), dir);
00581 }
00582
00583 template <typename To_Policy, typename From_Policy, typename Type>
00584 inline Result
00585 abs_float(Type& to, const Type from, Rounding_Dir) {
00586 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00587 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00588 to = std::abs(from);
00589 return V_EQ;
00590 }
00591
00592 template <typename To_Policy, typename From_Policy, typename Type>
00593 inline Result
00594 sqrt_float(Type& to, const Type from, Rounding_Dir dir) {
00595 if (To_Policy::check_nan_result && is_nan<From_Policy>(from))
00596 return assign_special<To_Policy>(to, VC_NAN, ROUND_IGNORE);
00597 if (To_Policy::check_sqrt_neg && from < 0)
00598 return assign_special<To_Policy>(to, V_SQRT_NEG, ROUND_IGNORE);
00599 prepare_inexact<To_Policy>(dir);
00600 if (fpu_direct_rounding(dir))
00601 to = std::sqrt(from);
00602 else {
00603 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00604 avoid_cse(from);
00605 to = std::sqrt(from);
00606 avoid_cse(to);
00607 fpu_restore_rounding_direction(old);
00608 }
00609 return result_relation<To_Policy>(dir);
00610 }
00611
00612 template <typename Policy, typename Type>
00613 inline Result
00614 sgn_float(const Type x) {
00615 return classify<Policy>(x, false, false, true);
00616 }
00617
00618 template <typename Policy1, typename Policy2, typename Type>
00619 inline Result
00620 cmp_float(const Type x, const Type y) {
00621 if (x > y)
00622 return V_GT;
00623 if (x < y)
00624 return V_LT;
00625 if (x == y)
00626 return V_EQ;
00627 return V_UNORD_COMP;
00628 }
00629
00630 template <typename To_Policy, typename From_Policy, typename To, typename From>
00631 inline Result
00632 assign_float_int_inexact(To& to, const From from, Rounding_Dir dir) {
00633 prepare_inexact<To_Policy>(dir);
00634 if (fpu_direct_rounding(dir))
00635 to = from;
00636 else {
00637 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00638 avoid_cse(from);
00639 to = from;
00640 avoid_cse(to);
00641 fpu_restore_rounding_direction(old);
00642 }
00643 return result_relation<To_Policy>(dir);
00644 }
00645
00646 template <typename To_Policy, typename From_Policy, typename To, typename From>
00647 inline Result
00648 assign_float_int(To& to, const From from, Rounding_Dir dir) {
00649 if (sizeof(From) * 8 > Float<To>::Binary::MANTISSA_BITS)
00650 return assign_float_int_inexact<To_Policy, From_Policy>(to, from, dir);
00651 else
00652 return assign_exact<To_Policy, From_Policy>(to, from, dir);
00653 }
00654
00655 template <typename Policy, typename T>
00656 inline Result
00657 set_neg_overflow_float(T& to, Rounding_Dir dir) {
00658 switch (round_dir(dir)) {
00659 case ROUND_UP:
00660 {
00661 Float<T> f;
00662 f.u.binary.set_max(true);
00663 to = f.value();
00664 return V_LT;
00665 }
00666 default:
00667 to = -HUGE_VAL;
00668 return V_GT;
00669 }
00670 }
00671
00672 template <typename Policy, typename T>
00673 inline Result
00674 set_pos_overflow_float(T& to, Rounding_Dir dir) {
00675 switch (round_dir(dir)) {
00676 case ROUND_DOWN:
00677 {
00678 Float<T> f;
00679 f.u.binary.set_max(false);
00680 to = f.value();
00681 return V_GT;
00682 }
00683 default:
00684 to = HUGE_VAL;
00685 return V_LT;
00686 }
00687 }
00688
00689 template <typename To_Policy, typename From_Policy, typename T>
00690 inline Result
00691 assign_float_mpz(T& to, const mpz_class& _from, Rounding_Dir dir)
00692 {
00693 mpz_srcptr from = _from.get_mpz_t();
00694 int sign = mpz_sgn(from);
00695 if (sign == 0) {
00696 to = 0;
00697 return V_EQ;
00698 }
00699 size_t exponent = mpz_sizeinbase(from, 2) - 1;
00700 if (exponent > static_cast<size_t>(Float<T>::Binary::EXPONENT_MAX)) {
00701 if (sign < 0)
00702 return set_neg_overflow_float<To_Policy>(to, dir);
00703 else
00704 return set_pos_overflow_float<To_Policy>(to, dir);
00705 }
00706 unsigned long zeroes = mpn_scan1(from->_mp_d, 0);
00707 size_t meaningful_bits = exponent - zeroes;
00708 mpz_t mantissa;
00709 mpz_init(mantissa);
00710 if (exponent > Float<T>::Binary::MANTISSA_BITS)
00711 mpz_tdiv_q_2exp(mantissa,
00712 from,
00713 exponent - Float<T>::Binary::MANTISSA_BITS);
00714 else
00715 mpz_mul_2exp(mantissa, from, Float<T>::Binary::MANTISSA_BITS - exponent);
00716 Float<T> f(to);
00717 f.u.binary.build(sign < 0, mantissa, exponent);
00718 mpz_clear(mantissa);
00719 to = f.value();
00720 if (meaningful_bits > Float<T>::Binary::MANTISSA_BITS) {
00721 if (sign < 0)
00722 return round_lt_float<To_Policy>(to, dir);
00723 else
00724 return round_gt_float<To_Policy>(to, dir);
00725 }
00726 return V_EQ;
00727 }
00728
00729 template <typename To_Policy, typename From_Policy, typename T>
00730 inline Result
00731 assign_float_mpq(T& to, const mpq_class& from, Rounding_Dir dir)
00732 {
00733 const mpz_class& _num = from.get_num();
00734 const mpz_class& _den = from.get_den();
00735 if (_den == 1)
00736 return assign_float_mpz<To_Policy, From_Policy>(to, _num, dir);
00737 mpz_srcptr num = _num.get_mpz_t();
00738 mpz_srcptr den = _den.get_mpz_t();
00739 int sign = mpz_sgn(num);
00740 signed long exponent = mpz_sizeinbase(num, 2) - mpz_sizeinbase(den, 2);
00741 if (exponent < Float<T>::Binary::EXPONENT_MIN_DENORM) {
00742 to = 0;
00743 inexact:
00744 if (sign < 0)
00745 return round_lt_float<To_Policy>(to, dir);
00746 else
00747 return round_gt_float<To_Policy>(to, dir);
00748 }
00749 if (exponent > static_cast<signed int>(Float<T>::Binary::EXPONENT_MAX + 1)) {
00750 overflow:
00751 if (sign < 0)
00752 return set_neg_overflow_float<To_Policy>(to, dir);
00753 else
00754 return set_pos_overflow_float<To_Policy>(to, dir);
00755 }
00756 unsigned int needed_bits = Float<T>::Binary::MANTISSA_BITS + 1;
00757 if (exponent < Float<T>::Binary::EXPONENT_MIN)
00758 needed_bits -= Float<T>::Binary::EXPONENT_MIN - exponent;
00759 mpz_t mantissa;
00760 mpz_init(mantissa);
00761 signed long shift = needed_bits - exponent;
00762 if (shift > 0) {
00763 mpz_mul_2exp(mantissa, num, shift);
00764 num = mantissa;
00765 }
00766 else if (shift < 0) {
00767 mpz_mul_2exp(mantissa, den, -shift);
00768 den = mantissa;
00769 }
00770 mpz_t r;
00771 mpz_init(r);
00772 mpz_tdiv_qr(mantissa, r, num, den);
00773 size_t bits = mpz_sizeinbase(mantissa, 2);
00774 bool inexact = (mpz_sgn(r) != 0);
00775 mpz_clear(r);
00776 if (bits == needed_bits + 1) {
00777 inexact = (inexact || mpz_odd_p(mantissa));
00778 mpz_div_2exp(mantissa, mantissa, 1);
00779 }
00780 else
00781 --exponent;
00782 if (exponent > static_cast<signed int>(Float<T>::Binary::EXPONENT_MAX)) {
00783 mpz_clear(mantissa);
00784 goto overflow;
00785 }
00786 else if (exponent < Float<T>::Binary::EXPONENT_MIN - 1) {
00787
00788 exponent = Float<T>::Binary::EXPONENT_MIN - 1;
00789 }
00790 Float<T> f(to);
00791 f.u.binary.build(sign < 0, mantissa, exponent);
00792 mpz_clear(mantissa);
00793 to = f.value();
00794 if (inexact)
00795 goto inexact;
00796 return V_EQ;
00797 }
00798
00799 template <typename To_Policy, typename From1_Policy, typename From2_Policy,
00800 typename Type>
00801 inline Result
00802 add_mul_float(Type& to, const Type x, const Type y, Rounding_Dir dir) {
00803 if (To_Policy::check_inf_mul_zero
00804 && ((x == 0 && is_inf_float<From2_Policy>(y)) ||
00805 (y == 0 && is_inf_float<From1_Policy>(x))))
00806 return assign_special<To_Policy>(to, V_INF_MUL_ZERO, ROUND_IGNORE);
00807
00808 prepare_inexact<To_Policy>(dir);
00809 if (fpu_direct_rounding(dir))
00810 to = fma(x, y, to);
00811 else if (fpu_inverse_rounding(dir))
00812 to = -limit_precision(fma(-x, y, -to));
00813 else {
00814 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00815 avoid_cse(x);
00816 avoid_cse(y);
00817 to = fma(x, y, to);
00818 avoid_cse(to);
00819 fpu_restore_rounding_direction(old);
00820 }
00821 if (To_Policy::check_nan_result && is_nan<To_Policy>(to))
00822 return VC_NAN;
00823 return result_relation<To_Policy>(dir);
00824 }
00825
00826 template <typename To_Policy, typename From1_Policy, typename From2_Policy, typename Type>
00827 inline Result
00828 sub_mul_float(Type& to, const Type x, const Type y, Rounding_Dir dir) {
00829 if (To_Policy::check_inf_mul_zero
00830 && ((x == 0 && is_inf_float<From2_Policy>(y)) ||
00831 (y == 0 && is_inf_float<From1_Policy>(x))))
00832 return assign_special<To_Policy>(to, V_INF_MUL_ZERO, ROUND_IGNORE);
00833
00834 prepare_inexact<To_Policy>(dir);
00835 if (fpu_direct_rounding(dir))
00836 to = fma(x, -y, to);
00837 else if (fpu_inverse_rounding(dir))
00838 to = -limit_precision(fma(x, y, -to));
00839 else {
00840 fpu_rounding_control_word_type old = fpu_save_rounding_direction(round_fpu_dir(dir));
00841 avoid_cse(x);
00842 avoid_cse(y);
00843 to = fma(x, -y, to);
00844 avoid_cse(to);
00845 fpu_restore_rounding_direction(old);
00846 }
00847 if (To_Policy::check_nan_result && is_nan<To_Policy>(to))
00848 return VC_NAN;
00849 return result_relation<To_Policy>(dir);
00850 }
00851
00852 template <typename Policy, typename Type>
00853 inline Result
00854 output_float(std::ostream& os, const Type from, const Numeric_Format&,
00855 Rounding_Dir) {
00856 if (from == 0)
00857 os << "0";
00858 else if (is_minf<Policy>(from))
00859 os << "-inf";
00860 else if (is_pinf<Policy>(from))
00861 os << "+inf";
00862 else if (is_nan<Policy>(from))
00863 os << "nan";
00864 else {
00865 int old_precision = os.precision(10000);
00866
00867
00868
00869
00870
00871 os << from;
00872 os.precision(old_precision);
00873 }
00874 return V_EQ;
00875 }
00876
00877 #if PPL_SUPPORTED_FLOAT
00878 SPECIALIZE_ASSIGN(assign_float_float_exact, float, float)
00879 #if PPL_SUPPORTED_DOUBLE
00880 SPECIALIZE_ASSIGN(assign_float_float, float, double)
00881 SPECIALIZE_ASSIGN(assign_float_float_exact, double, float)
00882 #endif
00883 #if PPL_SUPPORTED_LONG_DOUBLE
00884 SPECIALIZE_ASSIGN(assign_float_float, float, long double)
00885 SPECIALIZE_ASSIGN(assign_float_float_exact, long double, float)
00886 #endif
00887 #endif
00888
00889 #if PPL_SUPPORTED_DOUBLE
00890 SPECIALIZE_ASSIGN(assign_float_float_exact, double, double)
00891 #if PPL_SUPPORTED_LONG_DOUBLE
00892 SPECIALIZE_ASSIGN(assign_float_float, double, long double)
00893 SPECIALIZE_ASSIGN(assign_float_float_exact, long double, double)
00894 #endif
00895 #endif
00896
00897 #if PPL_SUPPORTED_LONG_DOUBLE
00898 SPECIALIZE_ASSIGN(assign_float_float_exact, long double, long double)
00899 #endif
00900
00901 #if PPL_SUPPORTED_FLOAT
00902 SPECIALIZE_CLASSIFY(classify_float, float)
00903 SPECIALIZE_IS_NAN(is_nan_float, float)
00904 SPECIALIZE_IS_MINF(is_minf_float, float)
00905 SPECIALIZE_IS_PINF(is_pinf_float, float)
00906 SPECIALIZE_ASSIGN_SPECIAL(assign_special_float, float)
00907 SPECIALIZE_ASSIGN(assign_float_int, float, signed char)
00908 SPECIALIZE_ASSIGN(assign_float_int, float, signed short)
00909 SPECIALIZE_ASSIGN(assign_float_int, float, signed int)
00910 SPECIALIZE_ASSIGN(assign_float_int, float, signed long)
00911 SPECIALIZE_ASSIGN(assign_float_int, float, signed long long)
00912 SPECIALIZE_ASSIGN(assign_float_int, float, unsigned char)
00913 SPECIALIZE_ASSIGN(assign_float_int, float, unsigned short)
00914 SPECIALIZE_ASSIGN(assign_float_int, float, unsigned int)
00915 SPECIALIZE_ASSIGN(assign_float_int, float, unsigned long)
00916 SPECIALIZE_ASSIGN(assign_float_int, float, unsigned long long)
00917 SPECIALIZE_ASSIGN(assign_float_mpz, float, mpz_class)
00918 SPECIALIZE_ASSIGN(assign_float_mpq, float, mpq_class)
00919 SPECIALIZE_COPY(copy_generic, float)
00920 SPECIALIZE_IS_INT(is_int_float, float)
00921 SPECIALIZE_FLOOR(floor_float, float, float)
00922 SPECIALIZE_CEIL(ceil_float, float, float)
00923 SPECIALIZE_TRUNC(trunc_float, float, float)
00924 SPECIALIZE_NEG(neg_float, float, float)
00925 SPECIALIZE_ABS(abs_float, float, float)
00926 SPECIALIZE_ADD(add_float, float, float, float)
00927 SPECIALIZE_SUB(sub_float, float, float, float)
00928 SPECIALIZE_MUL(mul_float, float, float, float)
00929 SPECIALIZE_DIV(div_float, float, float, float)
00930 SPECIALIZE_REM(rem_float, float, float, float)
00931 SPECIALIZE_MUL2EXP(mul2exp_float, float, float)
00932 SPECIALIZE_DIV2EXP(div2exp_float, float, float)
00933 SPECIALIZE_SQRT(sqrt_float, float, float)
00934 SPECIALIZE_GCD(gcd_exact, float, float, float)
00935 SPECIALIZE_GCDEXT(gcdext_exact, float, float, float, float, float)
00936 SPECIALIZE_LCM(lcm_gcd_exact, float, float, float)
00937 SPECIALIZE_SGN(sgn_float, float)
00938 SPECIALIZE_CMP(cmp_float, float, float)
00939 SPECIALIZE_ADD_MUL(add_mul_float, float, float, float)
00940 SPECIALIZE_SUB_MUL(sub_mul_float, float, float, float)
00941 SPECIALIZE_INPUT(input_generic, float)
00942 SPECIALIZE_OUTPUT(output_float, float)
00943 #endif
00944
00945 #if PPL_SUPPORTED_DOUBLE
00946 SPECIALIZE_CLASSIFY(classify_float, double)
00947 SPECIALIZE_IS_NAN(is_nan_float, double)
00948 SPECIALIZE_IS_MINF(is_minf_float, double)
00949 SPECIALIZE_IS_PINF(is_pinf_float, double)
00950 SPECIALIZE_ASSIGN_SPECIAL(assign_special_float, double)
00951 SPECIALIZE_ASSIGN(assign_float_int, double, signed char)
00952 SPECIALIZE_ASSIGN(assign_float_int, double, signed short)
00953 SPECIALIZE_ASSIGN(assign_float_int, double, signed int)
00954 SPECIALIZE_ASSIGN(assign_float_int, double, signed long)
00955 SPECIALIZE_ASSIGN(assign_float_int, double, signed long long)
00956 SPECIALIZE_ASSIGN(assign_float_int, double, unsigned char)
00957 SPECIALIZE_ASSIGN(assign_float_int, double, unsigned short)
00958 SPECIALIZE_ASSIGN(assign_float_int, double, unsigned int)
00959 SPECIALIZE_ASSIGN(assign_float_int, double, unsigned long)
00960 SPECIALIZE_ASSIGN(assign_float_int, double, unsigned long long)
00961 SPECIALIZE_ASSIGN(assign_float_mpz, double, mpz_class)
00962 SPECIALIZE_ASSIGN(assign_float_mpq, double, mpq_class)
00963 SPECIALIZE_COPY(copy_generic, double)
00964 SPECIALIZE_IS_INT(is_int_float, double)
00965 SPECIALIZE_FLOOR(floor_float, double, double)
00966 SPECIALIZE_CEIL(ceil_float, double, double)
00967 SPECIALIZE_TRUNC(trunc_float, double, double)
00968 SPECIALIZE_NEG(neg_float, double, double)
00969 SPECIALIZE_ABS(abs_float, double, double)
00970 SPECIALIZE_ADD(add_float, double, double, double)
00971 SPECIALIZE_SUB(sub_float, double, double, double)
00972 SPECIALIZE_MUL(mul_float, double, double, double)
00973 SPECIALIZE_DIV(div_float, double, double, double)
00974 SPECIALIZE_REM(rem_float, double, double, double)
00975 SPECIALIZE_MUL2EXP(mul2exp_float, double, double)
00976 SPECIALIZE_DIV2EXP(div2exp_float, double, double)
00977 SPECIALIZE_SQRT(sqrt_float, double, double)
00978 SPECIALIZE_GCD(gcd_exact, double, double, double)
00979 SPECIALIZE_GCDEXT(gcdext_exact, double, double, double, double, double)
00980 SPECIALIZE_LCM(lcm_gcd_exact, double, double, double)
00981 SPECIALIZE_SGN(sgn_float, double)
00982 SPECIALIZE_CMP(cmp_float, double, double)
00983 SPECIALIZE_ADD_MUL(add_mul_float, double, double, double)
00984 SPECIALIZE_SUB_MUL(sub_mul_float, double, double, double)
00985 SPECIALIZE_INPUT(input_generic, double)
00986 SPECIALIZE_OUTPUT(output_float, double)
00987 #endif
00988
00989 #if PPL_SUPPORTED_LONG_DOUBLE
00990 SPECIALIZE_CLASSIFY(classify_float, long double)
00991 SPECIALIZE_IS_NAN(is_nan_float, long double)
00992 SPECIALIZE_IS_MINF(is_minf_float, long double)
00993 SPECIALIZE_IS_PINF(is_pinf_float, long double)
00994 SPECIALIZE_ASSIGN_SPECIAL(assign_special_float, long double)
00995 SPECIALIZE_ASSIGN(assign_float_int, long double, signed char)
00996 SPECIALIZE_ASSIGN(assign_float_int, long double, signed short)
00997 SPECIALIZE_ASSIGN(assign_float_int, long double, signed int)
00998 SPECIALIZE_ASSIGN(assign_float_int, long double, signed long)
00999 SPECIALIZE_ASSIGN(assign_float_int, long double, signed long long)
01000 SPECIALIZE_ASSIGN(assign_float_int, long double, unsigned char)
01001 SPECIALIZE_ASSIGN(assign_float_int, long double, unsigned short)
01002 SPECIALIZE_ASSIGN(assign_float_int, long double, unsigned int)
01003 SPECIALIZE_ASSIGN(assign_float_int, long double, unsigned long)
01004 SPECIALIZE_ASSIGN(assign_float_int, long double, unsigned long long)
01005 SPECIALIZE_ASSIGN(assign_float_mpz, long double, mpz_class)
01006 SPECIALIZE_ASSIGN(assign_float_mpq, long double, mpq_class)
01007 SPECIALIZE_COPY(copy_generic, long double)
01008 SPECIALIZE_IS_INT(is_int_float, long double)
01009 SPECIALIZE_FLOOR(floor_float, long double, long double)
01010 SPECIALIZE_CEIL(ceil_float, long double, long double)
01011 SPECIALIZE_TRUNC(trunc_float, long double, long double)
01012 SPECIALIZE_NEG(neg_float, long double, long double)
01013 SPECIALIZE_ABS(abs_float, long double, long double)
01014 SPECIALIZE_ADD(add_float, long double, long double, long double)
01015 SPECIALIZE_SUB(sub_float, long double, long double, long double)
01016 SPECIALIZE_MUL(mul_float, long double, long double, long double)
01017 SPECIALIZE_DIV(div_float, long double, long double, long double)
01018 SPECIALIZE_REM(rem_float, long double, long double, long double)
01019 SPECIALIZE_MUL2EXP(mul2exp_float, long double, long double)
01020 SPECIALIZE_DIV2EXP(div2exp_float, long double, long double)
01021 SPECIALIZE_SQRT(sqrt_float, long double, long double)
01022 SPECIALIZE_GCD(gcd_exact, long double, long double, long double)
01023 SPECIALIZE_GCDEXT(gcdext_exact, long double, long double, long double,
01024 long double, long double)
01025 SPECIALIZE_LCM(lcm_gcd_exact, long double, long double, long double)
01026 SPECIALIZE_SGN(sgn_float, long double)
01027 SPECIALIZE_CMP(cmp_float, long double, long double)
01028 SPECIALIZE_ADD_MUL(add_mul_float, long double, long double, long double)
01029 SPECIALIZE_SUB_MUL(sub_mul_float, long double, long double, long double)
01030 SPECIALIZE_INPUT(input_generic, long double)
01031 SPECIALIZE_OUTPUT(output_float, long double)
01032 #endif
01033
01034 }
01035
01036 }
01037
01038 #endif // !defined(PPL_checked_int_inlines_hh)