Commit 30630a4d authored by Vincent Lefèvre's avatar Vincent Lefèvre
Browse files

[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)

[From SVN r11783 (branches/3.1)]
parent 4807f317
......@@ -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;
}
}
......
......@@ -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;
}
......@@ -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;
}
......
......@@ -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;
}
......@@ -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;
}
......
......@@ -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;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment