From 73f78618f790db13452a9be7e247dad47942816d Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 15 Dec 2008 15:11:46 +0000 Subject: The fix for bug 6604 "incorrect directed rounding in mpfr_strtofr" in r5664 was incomplete. Completed fix and added more testcases. [merged -r5665:5670 -r5749:5751 from trunk] git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/2.3@5752 280ebfd0-de03-0410-8827-d642c229c3f4 --- strtofr.c | 154 +++++++++++++++++++++++++++++++++++++------------------ tests/tstrtofr.c | 90 ++++++++++++++++++++++++++------ 2 files changed, 178 insertions(+), 66 deletions(-) diff --git a/strtofr.c b/strtofr.c index 9f22e10c2..076e15269 100644 --- a/strtofr.c +++ b/strtofr.c @@ -31,14 +31,24 @@ MA 02110-1301, USA. */ #define MPFR_MAX_BASE 62 struct parsed_string { - int negative; - int base; - unsigned char *mantissa, *mant; - size_t prec, alloc; - mp_exp_t exp_base, exp_bin; + int negative; /* non-zero iff the number is negative */ + int base; /* base of the string */ + unsigned char *mantissa; /* raw significand (without any point) */ + unsigned char *mant; /* stripped significand (without starting and + ending zeroes) */ + size_t prec; /* length of mant (zero for +/-0) */ + size_t alloc; /* allocation size of mantissa */ + mp_exp_t exp_base; /* number of digits before the point */ + mp_exp_t exp_bin; /* exponent in case base=2 or 16, and the pxxx + format is used (i.e., exponent is given in + base 10) */ }; -/* This table has been generated by the following program */ +/* This table has been generated by the following program. + For 2 <= b <= MPFR_MAX_BASE, + RedInvLog2Table[b-2][0] / RedInvLog2Table[b-2][1] + is an upper approximation of log(2)/log(b). +*/ static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = { {1UL, 1UL}, {53UL, 84UL}, @@ -441,7 +451,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) MPFR_ZIV_DECL (loop); MPFR_TMP_DECL (marker); - /* determine the minimal precision for the computation */ + /* initialize the working precision */ prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)); /* compute y as long as rounding is not possible */ @@ -449,60 +459,88 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) MPFR_ZIV_INIT (loop, prec); for (;;) { - /* initialize y to the value of 0.mant_s[0]...mant_s[pr-1] */ + /* Set y to the value of the ~prec most significant bits of pstr->mant + (as long as we guarantee correct rounding, we don't need to get + exactly prec bits). */ ysize = (prec - 1) / BITS_PER_MP_LIMB + 1; + /* prec bits corresponds to ysize limbs */ ysize_bits = ysize * BITS_PER_MP_LIMB; + /* and to ysize_bits >= prec > MPFR_PREC (x) bits */ y = (mp_limb_t*) MPFR_TMP_ALLOC ((2 * ysize + 1) * sizeof (mp_limb_t)); - y += ysize; - - /* required precision for str to have ~ysize limbs - We must have (2^(BITS_PER_MP_LIMB))^ysize ~= base^pstr_size - aka pstr_size = ceil (ysize*BITS_PER_MP_LIMB/log2(base)) - ysize ~ prec/BITS_PER_MP_LIMB and prec < Umax/2 => + y += ysize; /* y has (ysize+1) allocated limbs */ + + /* pstr_size is the number of characters we read in pstr->mant + to have at least ysize full limbs. + We must have base^(pstr_size-1) >= (2^(BITS_PER_MP_LIMB))^ysize + (in the worst case, the first digit is one and all others are zero). + i.e., pstr_size >= 1 + ysize*BITS_PER_MP_LIMB/log2(base) + Since ysize ~ prec/BITS_PER_MP_LIMB and prec < Umax/2 => ysize*BITS_PER_MP_LIMB can not overflow. - Compute pstr_size = ysize_bits * Num / Den - Where Num/Den ~ 1/log2(base) + We compute pstr_size = 1 + ceil(ysize_bits * Num / Den) + where Num/Den >= 1/log2(base) It is not exactly ceil(1/log2(base)) but could be one more (base 2) - Quite ugly since it tries to avoid overflow */ - pstr_size = ( ysize_bits / RedInvLog2Table[pstr->base-2][1] - * RedInvLog2Table[pstr->base-2][0] ) - + (( ysize_bits % RedInvLog2Table[pstr->base-2][1]) - * RedInvLog2Table[pstr->base-2][0] - / RedInvLog2Table[pstr->base-2][1] ) - + 1; - + Quite ugly since it tries to avoid overflow: + let Num = RedInvLog2Table[pstr->base-2][0] + and Den = RedInvLog2Table[pstr->base-2][1], + and ysize_bits = a*Den+b, + then ysize_bits * Num/Den = a*Num + (b * Num)/Den, + thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den + */ + { + unsigned long Num = RedInvLog2Table[pstr->base-2][0]; + unsigned long Den = RedInvLog2Table[pstr->base-2][1]; + pstr_size = ((ysize_bits / Den) * Num) + + (((ysize_bits % Den) * Num + Den - 1) / Den) + + 1; + } + + /* since pstr_size corresponds to at least ysize_bits full bits, + and ysize_bits > prec, the weight of the neglected part of + pstr->mant (if any) is < ulp(y) < ulp(x) */ + + /* if the number of wanted characters is more than what we have in + pstr->mant, round it down */ if (pstr_size >= pstr->prec) pstr_size = pstr->prec; MPFR_ASSERTD (pstr_size == (mp_exp_t) pstr_size); - /* convert str into binary */ + /* convert str into binary: note that pstr->mant is big endian, + thus no offset is needed */ real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base); - MPFR_ASSERTD ( real_ysize <= ysize+1); + MPFR_ASSERTD (real_ysize <= ysize+1); /* normalize y: warning we can get even get ysize+1 limbs! */ - MPFR_ASSERTD (y[real_ysize - 1] != 0); + MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */ count_leading_zeros (count, y[real_ysize - 1]); - exact = (real_ysize <= ysize); - if (exact != 0) /* shift y to the left in that case y should be exact */ + /* exact means that the number of limbs of the output of mpn_set_str + is less or equal to ysize */ + exact = real_ysize <= ysize; + if (exact) /* shift y to the left in that case y should be exact */ { + /* we have enough limbs to store {y, real_ysize} */ /* shift {y, num_limb} for count bits to the left */ if (count != 0) - mpn_lshift (y, y, real_ysize, count); - /* shift {y, num_limb} for (ysize-num_limb) limbs to the left */ + mpn_lshift (y + ysize - real_ysize, y, real_ysize, count); if (real_ysize != ysize) { - MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize); + if (count == 0) + MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize); MPN_ZERO (y, ysize - real_ysize); } /* for each bit shift decrease exponent of y */ /* (This should not overflow) */ exp = - ((ysize - real_ysize) * BITS_PER_MP_LIMB + count); } - else /* shift y for the right */ + else /* shift y to the right, by doing this we might lose some + bits from the result of mpn_set_str (in addition to the + characters neglected from pstr->mant) */ { /* shift {y, num_limb} for (BITS_PER_MP_LIMB - count) bits - to the right */ - mpn_rshift (y, y, real_ysize, BITS_PER_MP_LIMB - count); + to the right. FIXME: can we prove that count cannot be zero here, + since mpn_rshift does not accept a shift of BITS_PER_MP_LIMB? */ + MPFR_ASSERTD (count != 0); + exact = mpn_rshift (y, y, real_ysize, BITS_PER_MP_LIMB - count) == + MPFR_LIMB_ZERO; /* for each bit shift increase exponent of y */ exp = BITS_PER_MP_LIMB - count; } @@ -542,7 +580,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) result = y; err = 0; } - /* case pstr->exp_base > pstr_size */ + /* case non-power-of-two-base, and pstr->exp_base > pstr_size */ else if (pstr->exp_base > (mp_exp_t) pstr_size) { mp_limb_t *z; @@ -559,6 +597,13 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) goto overflow; exact = exact && (err == -1); + /* If exact is non zero, then z equals exactly the value of the + pstr_size most significant digits from pstr->mant, i.e., the + only difference can come from the neglected pstr->prec-pstr_size + least significant digits of pstr->mant. + If exact is zero, then z is rounded towards zero with respect + to that value. */ + /* multiply(y = 0.mant_s[0]...mant_s[pr-1])_base by base^(exp_s-g) */ mpn_mul_n (result, y, z, ysize); @@ -588,6 +633,8 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) exp --; } + /* if the low ysize limbs of {result, 2*ysize} are all zero, + then the result is still "exact" (if it was before) */ exact = exact && (mpn_scan1 (result, 0) >= (unsigned long) ysize_bits); result += ysize; @@ -598,7 +645,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) mp_limb_t *z; mp_exp_t exp_z; - result = (mp_limb_t*) MPFR_TMP_ALLOC ( (3*ysize+1) * BYTES_PER_MP_LIMB); + result = (mp_limb_t*) MPFR_TMP_ALLOC ((3*ysize+1) * BYTES_PER_MP_LIMB); /* set y to y * K^ysize */ y = y - ysize; /* we have allocated ysize limbs at y - ysize */ @@ -622,7 +669,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) /* compute y / z */ /* result will be put into result + n, and remainder into result */ mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y, - 2*ysize, z, ysize); + 2 * ysize, z, ysize); /* exp -= exp_z + ysize_bits with overflow checking and check that we can add/substract 2 to exp without overflow */ @@ -635,6 +682,8 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) MPFR_EXP_MIN+2, MPFR_EXP_MAX-2, goto overflow, goto underflow); err += 2; + /* if the remainder of the division is zero, then the result is + still "exact" if it was before */ exact = exact && (mpn_popcount (result, ysize) == 0); /* normalize result */ @@ -648,7 +697,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) } result += ysize; } - /* case exp_base = pstr_size */ + /* case exp_base = pstr_size: no multiplication or division needed */ else { /* base^(exp_s-pr) = 1 nothing to compute */ @@ -656,18 +705,16 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) err = 0; } - if (pstr_size < pstr->prec && exact - && ((rnd == GMP_RNDN) - || ((rnd == GMP_RNDD && pstr->negative) - || (rnd == GMP_RNDU && !pstr->negative)))) - { - /* Some significant digits might have been forgotten, if so result - is not exact. */ - size_t i; + /* If result is exact, we still have to consider the neglected part + of the input string. For a directed rounding, in that case we could + still correctly round, since the neglected part is less than + one ulp, but that would make the code more complex, and give a + speedup for rare cases only. */ + exact = exact && (pstr_size == pstr->prec); - for (i = pstr_size; exact && i < pstr->prec; i++) - exact = pstr->mant[i] == 0; - } + /* at this point, result is an approximation rounded towards zero + of the pstr_size most significant digits of pstr->mant, with + equality in case exact is non-zero. */ /* test if rounding is possible, and if so exit the loop */ if (exact || mpfr_can_round_raw (result, ysize, @@ -692,6 +739,13 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) exp ++; } + if (res == 0) /* fix ternary value */ + { + exact = exact && (pstr_size == pstr->prec); + if (!exact) + res = (pstr->negative) ? 1 : -1; + } + /* Set sign of x before exp since check_range needs a valid sign */ (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x); diff --git a/tests/tstrtofr.c b/tests/tstrtofr.c index 2243bef7e..48a0060f9 100644 --- a/tests/tstrtofr.c +++ b/tests/tstrtofr.c @@ -931,30 +931,88 @@ check_retval (void) } /* Bug found by Christoph Lauter (in mpfr_set_str). */ +static struct bug20081025_test { + mpfr_rnd_t rnd; + int inexact; + const char *str; + const char *binstr; +} Bug20081028Table[] = { + {GMP_RNDN, -1, "1.00000000000000000006", "1"}, + {GMP_RNDZ, -1, "1.00000000000000000006", "1"}, + {GMP_RNDU, +1, "1.00000000000000000006", + "10000000000000000000000000000001e-31"}, + {GMP_RNDD, -1, "1.00000000000000000006", "1"}, + + + {GMP_RNDN, +1, "-1.00000000000000000006", "-1"}, + {GMP_RNDZ, +1, "-1.00000000000000000006", "-1"}, + {GMP_RNDU, +1, "-1.00000000000000000006", "-1"}, + {GMP_RNDD, -1, "-1.00000000000000000006", + "-10000000000000000000000000000001e-31"}, + + {GMP_RNDN, +1, "0.999999999999999999999999999999999999999999999", "1"}, + {GMP_RNDZ, -1, "0.999999999999999999999999999999999999999999999", + "11111111111111111111111111111111e-32"}, + {GMP_RNDU, +1, "0.999999999999999999999999999999999999999999999", "1"}, + {GMP_RNDD, -1, "0.999999999999999999999999999999999999999999999", + "11111111111111111111111111111111e-32"}, + + {GMP_RNDN, -1, "-0.999999999999999999999999999999999999999999999", "-1"}, + {GMP_RNDZ, +1, "-0.999999999999999999999999999999999999999999999", + "-11111111111111111111111111111111e-32"}, + {GMP_RNDU, +1, "-0.999999999999999999999999999999999999999999999", + "-11111111111111111111111111111111e-32"}, + {GMP_RNDD, -1, "-0.999999999999999999999999999999999999999999999", "-1"} +}; + static void bug20081028 (void) { - mpfr_t x; - const char *s = "0.10000000000000000000000000000001E1"; - int res, err = 0; + int i; + int inexact, res; + mpfr_rnd_t rnd; + mpfr_t x, y; + char *s; mpfr_init2 (x, 32); - res = mpfr_strtofr (x, "1.00000000000000000006", NULL, 10, GMP_RNDU); - if (res <= 0) - { - printf ("Error in bug20081028: expected positive ternary value," - " got %d\n", res); - err = 1; - } - if (! mpfr_greater_p (x, __gmpfr_one)) + mpfr_init2 (y, 32); + for (i = 0 ; i < numberof (Bug20081028Table) ; i++) { - printf ("Error in bug20081028:\nExpected %s\nGot ", s); - mpfr_dump (x); - err = 1; + rnd = Bug20081028Table[i].rnd; + inexact = Bug20081028Table[i].inexact; + mpfr_set_str_binary (x, Bug20081028Table[i].binstr); + res = mpfr_strtofr (y, Bug20081028Table[i].str, &s, 10, rnd); + if (s == NULL || *s != 0) + { + printf ("Error in Bug20081028: strtofr didn't parse entire input\n" + "for (i=%d) Str=\"%s\"", i, Bug20081028Table[i].str); + exit (1); + } + if (! SAME_SIGN (res, inexact)) + { + printf ("Error in Bug20081028: expected %s ternary value, " + "got %d\nfor (i=%d) Rnd=%s Str=\"%s\"\n Set binary gives: ", + inexact > 0 ? "positive" : "negative", + res, i, mpfr_print_rnd_mode(rnd), Bug20081028Table[i].str); + mpfr_dump (x); + printf (" strtofr gives: "); + mpfr_dump (y); + exit (1); + } + if (mpfr_cmp (x, y)) + { + printf ("Error in Bug20081028: Results differ between strtofr and " + "set_binary\nfor (i=%d) Rnd=%s Str=\"%s\"\n" + " Set binary gives: ", + i, mpfr_print_rnd_mode(rnd), Bug20081028Table[i].str); + mpfr_dump (x); + printf (" strtofr gives: "); + mpfr_dump (y); + exit (1); + } } + mpfr_clear (y); mpfr_clear (x); - if (err) - exit (1); } int -- cgit v1.2.1