diff -Naur mpfr-2.1.1-p11/set_ld.c mpfr-2.1.1-p12/set_ld.c --- mpfr-2.1.1-p11/set_ld.c 2005-04-21 14:12:23.000000000 +0000 +++ mpfr-2.1.1-p12/set_ld.c 2005-07-29 14:22:46.000000000 +0000 @@ -46,7 +46,8 @@ mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode) { mpfr_t t, u; - int inexact, shift_exp = 0, inexact2 = 0; + int inexact, shift_exp = 0; + long double x; LONGDOUBLE_NAN_ACTION (d, goto nan); @@ -67,11 +68,14 @@ mpfr_init2 (t, MPFR_LDBL_MANT_DIG); mpfr_init2 (u, IEEE_DBL_MANT_DIG); - mpfr_set_ui (t, 0, GMP_RNDN); mpfr_save_emin_emax (); - while (d != (long double) 0.0) + + convert: + x = d; + mpfr_set_ui (t, 0, GMP_RNDN); + while (x != (long double) 0.0) { - if ((d > (long double) DBL_MAX) || ((-d) > (long double) DBL_MAX)) + if ((x > (long double) DBL_MAX) || ((-x) > (long double) DBL_MAX)) { long double div9, div10, div11, div12, div13; @@ -83,31 +87,31 @@ div11 = div10 * div10; /* 2^(2^11) */ div12 = div11 * div11; /* 2^(2^12) */ div13 = div12 * div12; /* 2^(2^13) */ - if (ABS(d) >= div13) + if (ABS(x) >= div13) { - d = d / div13; /* exact */ + x /= div13; /* exact */ shift_exp += 8192; } - if (ABS(d) >= div12) + if (ABS(x) >= div12) { - d = d / div12; /* exact */ + x /= div12; /* exact */ shift_exp += 4096; } - if (ABS(d) >= div11) + if (ABS(x) >= div11) { - d = d / div11; /* exact */ + x /= div11; /* exact */ shift_exp += 2048; } - if (ABS(d) >= div10) + if (ABS(x) >= div10) { - d = d / div10; /* exact */ + x /= div10; /* exact */ shift_exp += 1024; } - /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < d < 2^1024, + /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024, therefore we have one extra exponent reduction step */ - if (ABS(d) >= div9) + if (ABS(x) >= div9) { - d = d / div9; /* exact */ + x /= div9; /* exact */ shift_exp += 512; } } @@ -118,63 +122,85 @@ /* div9 = 2^(-2^9) */ div10 = div9 * div9; /* 2^(-2^10) */ div11 = div10 * div10; /* 2^(-2^11) if extended precision */ - /* since -DBL_MAX <= d <= DBL_MAX, the cast to double should not + /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not overflow here */ - inexact = mpfr_set_d (u, (double) d, GMP_RNDN); + inexact = mpfr_set_d (u, (double) x, GMP_RNDZ); MPFR_ASSERTD(inexact == 0); - if (d != (long double) 0.0 && - ABS(d) < div10 && + if (x != (long double) 0.0 && + ABS(x) < div10 && div11 != (long double) 0.0 && div11 / div10 == div10) /* possible underflow */ { long double div12, div13; - /* After the divisions, any bit of d must be >= div10, + /* After the divisions, any bit of x must be >= div10, hence the possible division by div9. */ div12 = div11 * div11; /* 2^(-2^12) */ div13 = div12 * div12; /* 2^(-2^13) */ - if (ABS(d) <= div13) + if (ABS(x) <= div13) { - d = d / div13; /* exact */ + x /= div13; /* exact */ shift_exp -= 8192; } - if (ABS(d) <= div12) + if (ABS(x) <= div12) { - d = d / div12; /* exact */ + x /= div12; /* exact */ shift_exp -= 4096; } - if (ABS(d) <= div11) + if (ABS(x) <= div11) { - d = d / div11; /* exact */ + x /= div11; /* exact */ shift_exp -= 2048; } - if (ABS(d) <= div10) + if (ABS(x) <= div10) { - d = d / div10; /* exact */ + x /= div10; /* exact */ shift_exp -= 1024; } - if (ABS(d) <= div9) + if (ABS(x) <= div9) { - d = d / div9; /* exact */ + x /= div9; /* exact */ shift_exp -= 512; } } else { - mpfr_add (t, t, u, GMP_RNDN); /* exact */ - if (!mpfr_number_p (t)) - break; - d = d - (long double) mpfr_get_d1 (u); /* exact */ + if (mpfr_add (t, t, u, GMP_RNDZ) != 0) + { + if (!mpfr_number_p (t)) + break; + /* Inexact. This cannot happen unless the C implementation + "lies" on the precision or when long doubles are + implemented with FP expansions like under Mac OS X. */ + if (MPFR_PREC (t) != MPFR_PREC (r) + 1) + { + /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX. + The precision MPFR_PREC (r) + 1 allows us to + deduce the rounding bit and the sticky bit. */ + mpfr_set_prec (t, MPFR_PREC (r) + 1); + goto convert; + } + else + { + mp_limb_t *tp; + int rb_mask; + + /* Since mpfr_add was inexact, the sticky bit is 1. */ + tp = MPFR_MANT (t); + rb_mask = MPFR_LIMB_ONE << + (BITS_PER_MP_LIMB - 1 - + (MPFR_PREC (r) & (BITS_PER_MP_LIMB - 1))); + if (rnd_mode == GMP_RNDN) + rnd_mode = (*tp & rb_mask) ^ MPFR_IS_NEG (t) ? + GMP_RNDU : GMP_RNDD; + *tp |= rb_mask; + break; + } + } + x -= (long double) mpfr_get_d1 (u); /* exact */ } } } - /* now t is exactly the input value d */ - inexact = mpfr_set (r, t, rnd_mode); - if (shift_exp > 0) - inexact2 = mpfr_mul_2exp (r, r, shift_exp, rnd_mode); - else if (shift_exp < 0) - inexact2 = mpfr_div_2exp (r, r, -shift_exp, rnd_mode); - if (inexact2) /* overflow */ - inexact = inexact2; + inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode); mpfr_clear (t); mpfr_clear (u); mpfr_restore_emin_emax (); diff -Naur mpfr-2.1.1-p11/tests/tset_ld.c mpfr-2.1.1-p12/tests/tset_ld.c --- mpfr-2.1.1-p11/tests/tset_ld.c 2005-04-21 13:00:38.000000000 +0000 +++ mpfr-2.1.1-p12/tests/tset_ld.c 2005-07-29 14:56:38.000000000 +0000 @@ -155,7 +155,7 @@ mpfr_test_init (); tests_machine_prec_long_double (); - mpfr_init2 (x, 113); + mpfr_init2 (x, MPFR_LDBL_MANT_DIG); /* check +0.0 and -0.0 */ d = 0.0; @@ -210,7 +210,7 @@ /* checks that 2^i, 2^i+1 and 2^i-1 are correctly converted */ d = 1.0; - for (i = 1; i <= 113; i++) + for (i = 1; i < MPFR_LDBL_MANT_DIG; i++) { d = 2.0 * d; /* d = 2^i */ check_set_get (d, x);