From c25dc97f895f808eeb2a45d363358c0f6760f339 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Thu, 12 Oct 2017 14:23:53 +0000 Subject: [src] Fixed the behavior of the mpfr_get_* functions in a very reduced exponent range. (merged changesets r11770,11774,11775,11777 and r11779 manually from the trunk) git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/3.1@11783 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/get_ld.c | 4 ++++ src/get_si.c | 15 +++++++++++++-- src/get_sj.c | 27 +++++++++++++++++---------- src/get_ui.c | 9 +++++++++ src/get_uj.c | 11 ++++++++++- src/get_z.c | 10 ++++++++++ 6 files changed, 63 insertions(+), 13 deletions(-) diff --git a/src/get_ld.c b/src/get_ld.c index b1d77c771..51bd01552 100644 --- a/src/get_ld.c +++ b/src/get_ld.c @@ -41,6 +41,9 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */ mpfr_t y, z; int sign; + MPFR_SAVE_EXPO_DECL (expo); + + MPFR_SAVE_EXPO_MARK (expo); /* first round x to the target long double precision, so that all subsequent operations are exact (this avoids double rounding @@ -103,6 +106,7 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode) } if (sign < 0) r = -r; + MPFR_SAVE_EXPO_FREE (expo); return r; } } diff --git a/src/get_si.c b/src/get_si.c index 134874c30..3afad5839 100644 --- a/src/get_si.c +++ b/src/get_si.c @@ -28,6 +28,7 @@ mpfr_get_si (mpfr_srcptr f, mpfr_rnd_t rnd) mpfr_prec_t prec; long s; mpfr_t x; + MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (!mpfr_fits_slong_p (f, rnd))) { @@ -39,14 +40,22 @@ mpfr_get_si (mpfr_srcptr f, mpfr_rnd_t rnd) if (MPFR_IS_ZERO (f)) return (long) 0; - /* determine prec of long */ - for (s = LONG_MIN, prec = 0; s != 0; s /= 2, prec++) + /* Determine the precision of long. |LONG_MIN| may have one more bit + as an integer, but in this case, this is a power of 2, thus fits + in a precision-prec floating-point number. */ + for (s = LONG_MAX, prec = 0; s != 0; s /= 2, prec++) { } + MPFR_SAVE_EXPO_MARK (expo); + /* first round to prec bits */ mpfr_init2 (x, prec); mpfr_rint (x, f, rnd); + /* The flags from mpfr_rint are the wanted ones. In particular, + it sets the inexact flag when necessary. */ + MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); + /* warning: if x=0, taking its exponent is illegal */ if (MPFR_UNLIKELY (MPFR_IS_ZERO(x))) s = 0; @@ -65,5 +74,7 @@ mpfr_get_si (mpfr_srcptr f, mpfr_rnd_t rnd) mpfr_clear (x); + MPFR_SAVE_EXPO_FREE (expo); + return s; } diff --git a/src/get_sj.c b/src/get_sj.c index 5e92199d0..f270f375c 100644 --- a/src/get_sj.c +++ b/src/get_sj.c @@ -35,6 +35,7 @@ mpfr_get_sj (mpfr_srcptr f, mpfr_rnd_t rnd) intmax_t r; mpfr_prec_t prec; mpfr_t x; + MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (!mpfr_fits_intmax_p (f, rnd))) { @@ -46,20 +47,24 @@ mpfr_get_sj (mpfr_srcptr f, mpfr_rnd_t rnd) if (MPFR_IS_ZERO (f)) return (intmax_t) 0; - /* determine the precision of intmax_t */ - for (r = MPFR_INTMAX_MIN, prec = 0; r != 0; r /= 2, prec++) + /* Determine the precision of intmax_t. |INTMAX_MIN| may have one + more bit as an integer, but in this case, this is a power of 2, + thus fits in a precision-prec floating-point number. */ + for (r = MPFR_INTMAX_MAX, prec = 0; r != 0; r /= 2, prec++) { } - /* Note: though INTMAX_MAX would have been sufficient for the conversion, - we chose INTMAX_MIN so that INTMAX_MIN - 1 is always representable in - precision prec; this is useful to detect overflows in MPFR_RNDZ (will - be needed later). */ - /* Now, r = 0. */ + MPFR_ASSERTD (r == 0); + + MPFR_SAVE_EXPO_MARK (expo); mpfr_init2 (x, prec); mpfr_rint (x, f, rnd); MPFR_ASSERTN (MPFR_IS_FP (x)); + /* The flags from mpfr_rint are the wanted ones. In particular, + it sets the inexact flag when necessary. */ + MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); + if (MPFR_NOTZERO (x)) { mp_limb_t *xp; @@ -67,15 +72,15 @@ mpfr_get_sj (mpfr_srcptr f, mpfr_rnd_t rnd) xp = MPFR_MANT (x); sh = MPFR_GET_EXP (x); - MPFR_ASSERTN ((mpfr_prec_t) sh <= prec); + MPFR_ASSERTN ((mpfr_prec_t) sh <= prec + 1); if (MPFR_INTMAX_MIN + MPFR_INTMAX_MAX != 0 - && MPFR_UNLIKELY ((mpfr_prec_t) sh == prec)) + && MPFR_UNLIKELY ((mpfr_prec_t) sh > prec)) { /* 2's complement and x <= INTMAX_MIN: in the case mp_limb_t has the same size as intmax_t, we cannot use the code in the for loop since the operations would be performed in unsigned arithmetic. */ - MPFR_ASSERTN (MPFR_IS_NEG (x) && (mpfr_powerof2_raw (x))); + MPFR_ASSERTN (MPFR_IS_NEG (x) && mpfr_powerof2_raw (x)); r = MPFR_INTMAX_MIN; } else if (MPFR_IS_POS (x)) @@ -117,6 +122,8 @@ mpfr_get_sj (mpfr_srcptr f, mpfr_rnd_t rnd) mpfr_clear (x); + MPFR_SAVE_EXPO_FREE (expo); + return r; } diff --git a/src/get_ui.c b/src/get_ui.c index 48a5f333f..3c084eaca 100644 --- a/src/get_ui.c +++ b/src/get_ui.c @@ -30,6 +30,7 @@ mpfr_get_ui (mpfr_srcptr f, mpfr_rnd_t rnd) mpfr_t x; mp_size_t n; mpfr_exp_t exp; + MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (!mpfr_fits_ulong_p (f, rnd))) { @@ -44,10 +45,16 @@ mpfr_get_ui (mpfr_srcptr f, mpfr_rnd_t rnd) for (s = ULONG_MAX, prec = 0; s != 0; s /= 2, prec ++) { } + MPFR_SAVE_EXPO_MARK (expo); + /* first round to prec bits */ mpfr_init2 (x, prec); mpfr_rint (x, f, rnd); + /* The flags from mpfr_rint are the wanted ones. In particular, + it sets the inexact flag when necessary. */ + MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); + /* warning: if x=0, taking its exponent is illegal */ if (MPFR_IS_ZERO(x)) s = 0; @@ -61,5 +68,7 @@ mpfr_get_ui (mpfr_srcptr f, mpfr_rnd_t rnd) mpfr_clear (x); + MPFR_SAVE_EXPO_FREE (expo); + return s; } diff --git a/src/get_uj.c b/src/get_uj.c index 805c8a1a8..8d5e11250 100644 --- a/src/get_uj.c +++ b/src/get_uj.c @@ -35,6 +35,7 @@ mpfr_get_uj (mpfr_srcptr f, mpfr_rnd_t rnd) uintmax_t r; mpfr_prec_t prec; mpfr_t x; + MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (!mpfr_fits_uintmax_p (f, rnd))) { @@ -50,12 +51,18 @@ mpfr_get_uj (mpfr_srcptr f, mpfr_rnd_t rnd) for (r = MPFR_UINTMAX_MAX, prec = 0; r != 0; r /= 2, prec++) { } - /* Now, r = 0. */ + MPFR_ASSERTD (r == 0); + + MPFR_SAVE_EXPO_MARK (expo); mpfr_init2 (x, prec); mpfr_rint (x, f, rnd); MPFR_ASSERTN (MPFR_IS_FP (x)); + /* The flags from mpfr_rint are the wanted ones. In particular, + it sets the inexact flag when necessary. */ + MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); + if (MPFR_NOTZERO (x)) { mp_limb_t *xp; @@ -76,6 +83,8 @@ mpfr_get_uj (mpfr_srcptr f, mpfr_rnd_t rnd) mpfr_clear (x); + MPFR_SAVE_EXPO_FREE (expo); + return r; } diff --git a/src/get_z.c b/src/get_z.c index 68249a465..3c4811735 100644 --- a/src/get_z.c +++ b/src/get_z.c @@ -29,6 +29,7 @@ mpfr_get_z (mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd) int inex; mpfr_t r; mpfr_exp_t exp; + MPFR_SAVE_EXPO_DECL (expo); if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (f))) { @@ -41,6 +42,8 @@ mpfr_get_z (mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd) return 0; } + MPFR_SAVE_EXPO_MARK (expo); + exp = MPFR_GET_EXP (f); /* if exp <= 0, then |f|<1, thus |o(f)|<=1 */ MPFR_ASSERTN (exp < 0 || exp <= MPFR_PREC_MAX); @@ -50,6 +53,11 @@ mpfr_get_z (mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd) MPFR_ASSERTN (inex != 1 && inex != -1); /* integral part of f is representable in r */ MPFR_ASSERTN (MPFR_IS_FP (r)); + + /* The flags from mpfr_rint are the wanted ones. In particular, + it sets the inexact flag when necessary. */ + MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); + exp = mpfr_get_z_2exp (z, r); if (exp >= 0) mpz_mul_2exp (z, z, exp); @@ -57,5 +65,7 @@ mpfr_get_z (mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd) mpz_fdiv_q_2exp (z, z, -exp); mpfr_clear (r); + MPFR_SAVE_EXPO_FREE (expo); + return inex; } -- cgit v1.2.1