Submitted By: Jim Gifford (jim at cross-lfs dot org) Date: 2009-01-05 Initial Package Version: 2.3.2 Origin: MPFR Website Upstream Status: Fixed Description: See http://www.mpfr.org Website Under Bugs diff -Naur mpfr-2.3.2.orig/strtofr.c mpfr-2.3.2/strtofr.c --- mpfr-2.3.2.orig/strtofr.c 2008-09-12 06:47:25.000000000 -0700 +++ mpfr-2.3.2/strtofr.c 2009-01-05 03:26:12.000000000 -0800 @@ -31,14 +31,24 @@ #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 @@ 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 @@ 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; + y += ysize; /* y has (ysize+1) allocated limbs */ - /* 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 => + /* 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 ((mp_exp_t) pstr_size == (mp_exp_t) pstr_size); + 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 @@ 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 @@ 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 @@ 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 @@ 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 @@ /* 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 @@ 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 @@ } 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,6 +705,17 @@ err = 0; } + /* 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); + + /* 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, (pstr->negative) ? -1 : 1, @@ -679,6 +739,13 @@ 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 -Naur mpfr-2.3.2.orig/tests/tset_str.c mpfr-2.3.2/tests/tset_str.c --- mpfr-2.3.2.orig/tests/tset_str.c 2008-09-12 06:47:25.000000000 -0700 +++ mpfr-2.3.2/tests/tset_str.c 2009-01-05 03:25:12.000000000 -0800 @@ -81,6 +81,24 @@ mpfr_clear (a); } +/* Bug found by Christoph Lauter. */ +static void +bug20081028 (void) +{ + mpfr_t x; + const char *s = "0.10000000000000000000000000000001E1"; + + mpfr_init2 (x, 32); + mpfr_set_str (x, "1.00000000000000000006", 10, GMP_RNDU); + if (! mpfr_greater_p (x, __gmpfr_one)) + { + printf ("Error in bug20081028:\nExpected %s\nGot ", s); + mpfr_dump (x); + exit (1); + } + mpfr_clear (x); +} + int main (int argc, char *argv[]) { @@ -845,6 +863,7 @@ mpfr_clear (y); check_underflow (); + bug20081028 (); tests_end_mpfr (); return 0; diff -Naur mpfr-2.3.2.orig/tests/tstrtofr.c mpfr-2.3.2/tests/tstrtofr.c --- mpfr-2.3.2.orig/tests/tstrtofr.c 2008-09-12 06:47:25.000000000 -0700 +++ mpfr-2.3.2/tests/tstrtofr.c 2009-01-05 03:26:12.000000000 -0800 @@ -27,28 +27,7 @@ #include "mpfr-test.h" -static void check_reftable (void); -static void check_special (void); -static void check_retval (void); -static void check_overflow (void); -static void check_parse (void); - -int -main (int argc, char *argv[]) -{ - tests_start_mpfr (); - - check_special(); - check_reftable (); - check_parse (); - check_overflow (); - check_retval (); - - tests_end_mpfr (); - return 0; -} - -void +static void check_special (void) { mpfr_t x, y; @@ -551,8 +530,7 @@ "1.001000010110011011000101100000101111101001100011101101001111110111000010010110010001100e-16920"} }; - -void +static void check_reftable () { int i, base; @@ -597,7 +575,7 @@ mpfr_clear (x); } -void +static void check_parse (void) { mpfr_t x; @@ -951,3 +929,104 @@ mpfr_clear (x); } + +/* 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) +{ + int i; + int inexact, res; + mpfr_rnd_t rnd; + mpfr_t x, y; + char *s; + + mpfr_init2 (x, 32); + mpfr_init2 (y, 32); + for (i = 0 ; i < numberof (Bug20081028Table) ; i++) + { + 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); +} + +int +main (int argc, char *argv[]) +{ + tests_start_mpfr (); + + check_special(); + check_reftable (); + check_parse (); + check_overflow (); + check_retval (); + bug20081028 (); + + tests_end_mpfr (); + return 0; +}