diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2019-01-31 12:05:55.995705752 +0000 +++ mpfr-4.0.1-b/PATCHES 2019-01-31 12:05:56.091705514 +0000 @@ -0,0 +1 @@ +strtofr diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-12-31 11:18:05.621387341 +0000 +++ mpfr-4.0.1-b/VERSION 2019-01-31 12:05:56.091705514 +0000 @@ -1 +1 @@ -4.0.1-p14 +4.0.1-p15 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-12-31 11:18:05.617387407 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2019-01-31 12:05:56.079705544 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p14" +#define MPFR_VERSION_STRING "4.0.1-p15" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/strtofr.c mpfr-4.0.1-b/src/strtofr.c --- mpfr-4.0.1-a/src/strtofr.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/strtofr.c 2019-01-31 12:05:56.051705613 +0000 @@ -36,79 +36,80 @@ allocated for the mantissa field. */ size_t prec; /* length of mant (zero for +/-0) */ size_t alloc; /* allocation size of mantissa */ - mpfr_exp_t exp_base; /* number of digits before the point */ - mpfr_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) */ + mpfr_exp_t exp_base; /* number of digits before the point, + exponent + except in case of binary exponent (exp_bin) */ + mpfr_exp_t exp_bin; /* binary exponent of the pxxx format for + base = 2 or 16 */ }; /* 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). + is an upper approximation to log(2)/log(b), no larger than 1. + Note: these numbers must fit on 16 bits, thus unsigned int is OK. */ -static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = { - {1UL, 1UL}, - {53UL, 84UL}, - {1UL, 2UL}, - {4004UL, 9297UL}, - {53UL, 137UL}, - {2393UL, 6718UL}, - {1UL, 3UL}, - {665UL, 2108UL}, - {4004UL, 13301UL}, - {949UL, 3283UL}, - {53UL, 190UL}, - {5231UL, 19357UL}, - {2393UL, 9111UL}, - {247UL, 965UL}, - {1UL, 4UL}, - {4036UL, 16497UL}, - {665UL, 2773UL}, - {5187UL, 22034UL}, - {4004UL, 17305UL}, - {51UL, 224UL}, - {949UL, 4232UL}, - {3077UL, 13919UL}, - {53UL, 243UL}, - {73UL, 339UL}, - {5231UL, 24588UL}, - {665UL, 3162UL}, - {2393UL, 11504UL}, - {4943UL, 24013UL}, - {247UL, 1212UL}, - {3515UL, 17414UL}, - {1UL, 5UL}, - {4415UL, 22271UL}, - {4036UL, 20533UL}, - {263UL, 1349UL}, - {665UL, 3438UL}, - {1079UL, 5621UL}, - {5187UL, 27221UL}, - {2288UL, 12093UL}, - {4004UL, 21309UL}, - {179UL, 959UL}, - {51UL, 275UL}, - {495UL, 2686UL}, - {949UL, 5181UL}, - {3621UL, 19886UL}, - {3077UL, 16996UL}, - {229UL, 1272UL}, - {53UL, 296UL}, - {109UL, 612UL}, - {73UL, 412UL}, - {1505UL, 8537UL}, - {5231UL, 29819UL}, - {283UL, 1621UL}, - {665UL, 3827UL}, - {32UL, 185UL}, - {2393UL, 13897UL}, - {1879UL, 10960UL}, - {4943UL, 28956UL}, - {409UL, 2406UL}, - {247UL, 1459UL}, - {231UL, 1370UL}, - {3515UL, 20929UL} }; +static const unsigned int RedInvLog2Table[MPFR_MAX_BASE-1][2] = { + {1, 1}, + {53, 84}, + {1, 2}, + {4004, 9297}, + {53, 137}, + {2393, 6718}, + {1, 3}, + {665, 2108}, + {4004, 13301}, + {949, 3283}, + {53, 190}, + {5231, 19357}, + {2393, 9111}, + {247, 965}, + {1, 4}, + {4036, 16497}, + {665, 2773}, + {5187, 22034}, + {4004, 17305}, + {51, 224}, + {949, 4232}, + {3077, 13919}, + {53, 243}, + {73, 339}, + {5231, 24588}, + {665, 3162}, + {2393, 11504}, + {4943, 24013}, + {247, 1212}, + {3515, 17414}, + {1, 5}, + {4415, 22271}, + {4036, 20533}, + {263, 1349}, + {665, 3438}, + {1079, 5621}, + {5187, 27221}, + {2288, 12093}, + {4004, 21309}, + {179, 959}, + {51, 275}, + {495, 2686}, + {949, 5181}, + {3621, 19886}, + {3077, 16996}, + {229, 1272}, + {53, 296}, + {109, 612}, + {73, 412}, + {1505, 8537}, + {5231, 29819}, + {283, 1621}, + {665, 3827}, + {32, 185}, + {2393, 13897}, + {1879, 10960}, + {4943, 28956}, + {409, 2406}, + {247, 1459}, + {231, 1370}, + {3515, 20929} }; #if 0 #define N 8 int main () @@ -377,6 +378,13 @@ res = 1; MPFR_ASSERTD (pstr->exp_base >= 0); + /* FIXME: In the code below (both cases), if the exponent from the + string is large, it will be replaced by MPFR_EXP_MIN or MPFR_EXP_MAX, + i.e. it will have a different value. This may not change the result + in most cases, but there is no guarantee on very long strings when + mpfr_exp_t is a 32-bit type, as the exponent could be brought back + to the current exponent range. */ + /* an optional exponent (e or E, p or P, @) */ if ( (*str == '@' || (base <= 10 && (*str == 'e' || *str == 'E'))) && (!isspace((unsigned char) str[1])) ) @@ -445,112 +453,174 @@ static int parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd) { - mpfr_prec_t prec; + mpfr_prec_t precx, prec, ysize_bits, pstr_size; mpfr_exp_t exp; - mpfr_exp_t ysize_bits; - mp_limb_t *y, *result; + mp_limb_t *result; int count, exact; - size_t pstr_size; - mp_size_t ysize, real_ysize; + mp_size_t ysize, real_ysize, diff_ysize; int res, err; + const int extra_limbs = GMP_NUMB_BITS >= 12 ? 1 : 2; /* see below */ MPFR_ZIV_DECL (loop); MPFR_TMP_DECL (marker); /* initialize the working precision */ - prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)); + precx = MPFR_GET_PREC (x); + prec = precx + MPFR_INT_CEIL_LOG2 (precx); - /* compute the value y of the leading characters as long as rounding is not - possible */ + /* Compute the value y of the leading characters as long as rounding is not + possible. + Note: We have some integer overflow checking using MPFR_EXP_MIN and + MPFR_EXP_MAX in this loop. Thanks to the large margin between these + extremal values of the mpfr_exp_t type and the valid minimum/maximum + exponents, such integer overflows would correspond to real underflow + or overflow on the result (possibly except in huge precisions, which + are disregarded here; anyway, in practice, such issues could occur + only with 32-bit precision and exponent types). Such checks could be + extended to real early underflow/overflow checking, in order to avoid + useless computations in such cases; in such a case, be careful that + the approximation errors need to be taken into account. */ MPFR_TMP_MARK(marker); MPFR_ZIV_INIT (loop, prec); for (;;) { - /* 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). */ + mp_limb_t *y0, *y; + + /* y will be regarded as a number with precision prec. */ ysize = MPFR_PREC2LIMBS (prec); /* prec bits corresponds to ysize limbs */ - ysize_bits = ysize * GMP_NUMB_BITS; - /* and to ysize_bits >= prec > MPFR_PREC (x) bits */ - /* we need to allocate one more limb to work around bug + ysize_bits = (mpfr_prec_t) ysize * GMP_NUMB_BITS; + MPFR_ASSERTD (ysize_bits >= prec); + /* and to ysize_bits >= prec > precx bits. */ + /* We need to allocate one more limb as specified by mpn_set_str + (a limb may be written in rp[rn]). Note that the manual of GMP + up to 5.1.3 was incorrect on this point. + See the following discussion: https://gmplib.org/list-archives/gmp-bugs/2013-December/003267.html */ - y = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 2); - y += ysize; /* y has (ysize+2) allocated limbs */ + y0 = MPFR_TMP_LIMBS_ALLOC (2 * ysize + extra_limbs + 1); + y = y0 + ysize; /* y has (ysize + extra_limbs + 1) allocated limbs */ - /* pstr_size is the number of characters we read in pstr->mant - to have at least ysize full limbs. + /* pstr_size is the number of bytes we want to read from pstr->mant + to fill at least ysize full limbs with mpn_set_str. We must have base^(pstr_size-1) >= (2^(GMP_NUMB_BITS))^ysize (in the worst case, the first digit is one and all others are zero). i.e., pstr_size >= 1 + ysize*GMP_NUMB_BITS/log2(base) Since ysize ~ prec/GMP_NUMB_BITS and prec < Umax/2 => ysize*GMP_NUMB_BITS can not overflow. 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) + where 1/log2(base) <= Num/Den <= 1 + It is not exactly ceil(1/log2(base)) but could be one more (base 2). 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 + + Note: denoting m = pstr_size and n = ysize_bits, assuming we have + m = 1 + ceil(n/log2(b)), i.e., b^(m-1) >= 2^n > b^(m-2), then + b^(m-1)/2^n < b, and since we consider m characters of the input, + the corresponding part is less than b^m < b^2*2^n. + This implies that if b^2 < 2^GMP_NUMB_BITS, which for b <= 62 holds + for GMP_NUMB_BITS >= 12, we have real_ysize <= ysize+1 below + (this also implies that for GMP_NUMB_BITS >= 13, the number of bits + of y[real_ysize-1] below is less than GMP_NUMB_BITS, thus + count < GMP_NUMB_BITS). + Warning: for GMP_NUMB_BITS=8, we can have real_ysize = ysize + 2! + Hence the allocation above for ysize + extra_limbs limbs. */ { - 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) + unsigned int Num = RedInvLog2Table[pstr->base-2][0]; + unsigned int Den = RedInvLog2Table[pstr->base-2][1]; + MPFR_ASSERTD (Num <= Den && Den <= 65535); /* thus no overflow */ + pstr_size = (ysize_bits / Den) * Num + + ((unsigned long) (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) */ + /* Since pstr_size corresponds to at least ysize_bits 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) + /* If the number of wanted bytes is more than what is available + in pstr->mant, i.e. pstr->prec, reduce it to pstr->prec. */ + if (pstr_size > pstr->prec) pstr_size = pstr->prec; - MPFR_ASSERTD (pstr_size == (mpfr_exp_t) pstr_size); - /* convert str into binary: note that pstr->mant is big endian, - thus no offset is needed */ + /* Convert str (potentially truncated to pstr_size) 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); - /* normalize y: warning we can even get ysize+1 limbs! */ + /* See above for the explanation of the following assertion. */ + MPFR_ASSERTD (real_ysize <= ysize + extra_limbs); + + /* The Boolean "exact" will attempt to track exactness of the result: + If it is true, then this means that the result is exact, allowing + termination, even though the rounding test may not succeed. + Conversely, if the result is exact, then "exact" will not + necessarily be true at the end of the Ziv loop, but we will need + to make sure that at some point, "exact" will be true in order to + guarantee termination. FIXME: check that. */ + /* First, consider the part of the input string that has been ignored. + Note that the trailing zeros have been removed in parse_string, so + that if something has been ignored, it must be non-zero. */ + exact = pstr_size == pstr->prec; + + /* Normalize y and set the initial value of its exponent exp, which + is 0 when y is not shifted. + Since pstr->mant was normalized, mpn_set_str guarantees that + the most significant limb is non-zero. */ MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */ count_leading_zeros (count, y[real_ysize - 1]); - /* 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 */ + diff_ysize = ysize - real_ysize; + MPFR_LOG_MSG (("diff_ysize = %ld\n", (long) diff_ysize)); + if (diff_ysize >= 0) { - /* we have enough limbs to store {y, real_ysize} */ - /* shift {y, num_limb} for count bits to the left */ + /* We have enough limbs to store {y, real_ysize} exactly + in {y, ysize}, so that we can do a left shift, without + losing any information ("exact" will not change). */ if (count != 0) - mpn_lshift (y + ysize - real_ysize, y, real_ysize, count); - if (real_ysize != ysize) + mpn_lshift (y + diff_ysize, y, real_ysize, count); + if (diff_ysize > 0) { if (count == 0) - mpn_copyd (y + ysize - real_ysize, y, real_ysize); - MPN_ZERO (y, ysize - real_ysize); + mpn_copyd (y + diff_ysize, y, real_ysize); + MPN_ZERO (y, diff_ysize); } - /* for each bit shift decrease exponent of y */ - /* (This should not overflow) */ - exp = - ((ysize - real_ysize) * GMP_NUMB_BITS + count); + /* exp = negation of the total shift count, avoiding overflows. */ + exp = - ((mpfr_exp_t) diff_ysize * GMP_NUMB_BITS + count); } - 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) */ + else { - /* shift {y, num_limb} for (GMP_NUMB_BITS - count) bits - to the right. FIXME: can we prove that count cannot be zero here, - since mpn_rshift does not accept a shift of GMP_NUMB_BITS? */ - MPFR_ASSERTD (count != 0); - exact = mpn_rshift (y, y, real_ysize, GMP_NUMB_BITS - count) == - MPFR_LIMB_ZERO; - /* for each bit shift increase exponent of y */ - exp = GMP_NUMB_BITS - count; + /* Shift {y, real_ysize} for (GMP_NUMB_BITS - count) bits to the + right, and put the ysize most significant limbs into {y, ysize}. + We have either real_ysize = ysize + 1 or real_ysize = ysize + 2 + (only possible with extra_limbs == 2). */ + MPFR_ASSERTD (diff_ysize == -1 || + (extra_limbs == 2 && diff_ysize == -2)); + if (count != 0) + { + /* Before doing the shift, consider the limb that will entirely + be lost if real_ysize = ysize + 2. */ + exact = exact && (diff_ysize == -1 || y[0] == MPFR_LIMB_ZERO); + /* mpn_rshift allows overlap, provided destination <= source */ + /* FIXME: The bits lost due to mpn_rshift are not taken + into account in the error analysis below! */ + if (mpn_rshift (y, y - (diff_ysize + 1), real_ysize, + GMP_NUMB_BITS - count) != MPFR_LIMB_ZERO) + exact = 0; /* some non-zero bits have been shifted out */ + } + else + { + /* the case real_ysize = ysize + 2 with count = 0 cannot happen + even with GMP_NUMB_BITS = 8 since 62^2 < 256^2/2 */ + MPFR_ASSERTD (diff_ysize == -1); + exact = exact && y[0] == MPFR_LIMB_ZERO; + /* copy {y+real_ysize-ysize, ysize} to {y, ysize} */ + mpn_copyi (y, y + 1, real_ysize - 1); + } + /* exp = shift count */ + /* TODO: add some explanations about what exp means exactly. */ + exp = GMP_NUMB_BITS * (- diff_ysize) - count; } /* compute base^(exp_base - pstr_size) on n limbs */ @@ -560,6 +630,8 @@ int pow2; mpfr_exp_t tmp; + MPFR_LOG_MSG (("case 1 (base = power of 2)\n", 0)); + count_leading_zeros (pow2, (mp_limb_t) pstr->base); pow2 = GMP_NUMB_BITS - pow2 - 1; /* base = 2^pow2 */ MPFR_ASSERTD (0 < pow2 && pow2 <= 5); @@ -596,10 +668,12 @@ mp_limb_t *z; mpfr_exp_t exp_z; + MPFR_LOG_MSG (("case 2 (exp_base > pstr_size)\n", 0)); + result = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1); /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */ - z = y - ysize; + z = y0; /* NOTE: exp_base-pstr_size can't overflow since pstr_size > 0 */ err = mpfr_mpn_exp (z, &exp_z, pstr->base, pstr->exp_base - pstr_size, ysize); @@ -646,8 +720,7 @@ /* 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); + exact = exact && (mpn_scan1 (result, 0) >= ysize_bits); result += ysize; } /* case exp_base < pstr_size */ @@ -656,11 +729,12 @@ mp_limb_t *z; mpfr_exp_t exp_z; + MPFR_LOG_MSG (("case 3 (exp_base < pstr_size)\n", 0)); + result = MPFR_TMP_LIMBS_ALLOC (3 * ysize + 1); - /* set y to y * K^ysize */ - y = y - ysize; /* we have allocated ysize limbs at y - ysize */ - MPN_ZERO (y, ysize); + /* y0 = y * K^ysize */ + MPN_ZERO (y0, ysize); /* pstr_size - pstr->exp_base can overflow */ MPFR_SADD_OVERFLOW (exp_z, (mpfr_exp_t) pstr_size, -pstr->exp_base, @@ -668,33 +742,35 @@ MPFR_EXP_MIN, MPFR_EXP_MAX, goto underflow, goto overflow); - /* (z, exp_z) = base^(exp_base-pstr_size) */ + /* (z, exp_z) = base^(pstr_size - exp_base) */ z = result + 2*ysize + 1; err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize); - /* Since we want y/z rounded toward zero, we must get an upper - bound of z. If err >= 0, the error on z is bounded by 2^err. */ - if (err >= 0) - { - mp_limb_t cy; - unsigned long h = err / GMP_NUMB_BITS; - unsigned long l = err - h * GMP_NUMB_BITS; - if (h >= ysize) /* not enough precision in z */ - goto next_loop; - cy = mpn_add_1 (z, z, ysize - h, MPFR_LIMB_ONE << l); - if (cy != 0) /* the code below requires z on ysize limbs */ - goto next_loop; - } - exact = exact && (err == -1); + /* Now {z, ysize} * 2^(exp_z_out - ysize_bits) is an approximation + to base^exp_z_in (denoted b^e below), rounded toward zero, with: + * if err = -1, the result is exact; + * if err = -2, an overflow occurred in the computation of exp_z; + * otherwise the error is bounded by 2^err ulps. + Thus the exact value of b^e is between z and z + 2^err, where + z is {z, ysize} properly scaled by a power of 2. Then the error + will be: + y/b^e - trunc(y/z) = eps1 + eps2 + with + eps1 = y/b^e - y/z <= 0 + eps2 = y/z - trunc(y/z) >= 0 + thus the errors will (partly) compensate, giving a bound + max(|eps1|,|eps2|). + In addition, there is a 3rd error eps3 since y might be the + conversion of only a part of the character string, and/or y + might be truncated by the mpn_rshift call above: + eps3 = exact_y/b^e - y/b^e >= 0. + */ if (err == -2) goto underflow; /* FIXME: Sure? */ - if (err == -1) - err = 0; - - /* 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); + else if (err == -1) + err = 0; /* see the note below */ + else + exact = 0; /* exp -= exp_z + ysize_bits with overflow checking and check that we can add/subtract 2 to exp without overflow */ @@ -706,7 +782,59 @@ mpfr_exp_t, mpfr_uexp_t, MPFR_EXP_MIN+2, MPFR_EXP_MAX-2, goto overflow, goto underflow); - err += 2; + + /* Compute the integer division y/z rounded toward zero. + The quotient will be put at result + ysize (size: ysize + 1), + and the remainder at result (size: ysize). + Both the dividend {y, 2*ysize} and the divisor {z, ysize} are + normalized, i.e., the most significant bit of their most + significant limb is 1. */ + MPFR_ASSERTD (MPFR_LIMB_MSB (y0[2 * ysize - 1]) != 0); + MPFR_ASSERTD (MPFR_LIMB_MSB (z[ysize - 1]) != 0); + mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y0, + 2 * ysize, z, ysize); + + /* The truncation error of the mpn_tdiv_qr call (eps2 above) is at + most 1 ulp. Idem for the error eps3, which has the same sign, + thus eps2 + eps3 <= 2 ulps. + FIXME: For eps3, this is not obvious and should be explained. + For the error eps1 coming from the approximation to b^e, + we have (still up to a power-of-2 normalization): + y/z - y/b^e = y * (b^e-z) / (z * b^e) <= y * 2^err / (z * b^e). + We have to convert that error in terms of ulp(trunc(y/z)). + We first have ulp(trunc(y/z)) = ulp(y/z). + + FIXME: There must be some discussion about the exponents, + because up to a power of 2, 1/2 <= |y/z| < 1 and + 1 <= |y/z| < 2 are equivalent and give no information. + Moreover 1/2 <= b^e < 1 has not been explained and may + hide mistakes since one may have 1/2 <= z < 1 < b^e. + + Since both y and z are normalized, the quotient + {result+ysize, ysize+1} has exactly ysize limbs, plus maybe one + bit (this corresponds to the MPFR_ASSERTD below): + * if the quotient has exactly ysize limbs, then 1/2 <= |y/z| < 1 + (up to a power of 2) and since 1/2 <= b^e < 1, the error is at + most 2^(err+1) ulps; + * if the quotient has one extra bit, then 1 <= |y/z| < 2 + (up to a power of 2) and since 1/2 <= b^e < 1, the error is at + most 2^(err+2) ulps; but since we will shift the result right + below by one bit, the final error will be at most 2^(err+1) ulps + too. + + Thus the error is: + * at most 2^(err+1) ulps for eps1 + * at most 2 ulps for eps2 + eps3, which is of opposite sign + and we can bound the error by 2^(err+1) ulps in all cases. + + Note: If eps1 was 0, the error would be bounded by 2 ulps, + thus replacing err = -1 by err = 0 above was the right thing + to do, since 2^(0+1) = 2. + */ + MPFR_ASSERTD (result[2 * ysize] <= 1); + + err += 1; /* see above for the explanation of the +1 term */ + /* 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); @@ -726,17 +854,15 @@ /* case exp_base = pstr_size: no multiplication or division needed */ else { + MPFR_LOG_MSG (("case 4 (exp_base = pstr_size)\n", 0)); + /* base^(exp-pr) = 1 nothing to compute */ result = y; 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); + MPFR_LOG_MSG (("exact = %d, err = %d, precx = %Pu\n", + exact, err, precx)); /* at this point, result is an approximation rounded toward zero of the pstr_size most significant digits of pstr->mant, with @@ -744,24 +870,22 @@ /* test if rounding is possible, and if so exit the loop. Note: we also need to be able to determine the correct ternary value, - thus we use the MPFR_PREC(x) + (rnd == MPFR_RNDN) trick. + thus we use the precx + (rnd == MPFR_RNDN) trick. For example if result = xxx...xxx111...111 and rnd = RNDN, then we know the correct rounding is xxx...xx(x+1), but we cannot know the correct ternary value. */ if (exact || mpfr_round_p (result, ysize, ysize_bits - err - 1, - MPFR_PREC(x) + (rnd == MPFR_RNDN))) + precx + (rnd == MPFR_RNDN))) break; - next_loop: /* update the prec for next loop */ MPFR_ZIV_NEXT (loop, prec); } /* loop */ MPFR_ZIV_FREE (loop); /* round y */ - if (mpfr_round_raw (MPFR_MANT (x), result, - ysize_bits, - pstr->negative, MPFR_PREC(x), rnd, &res )) + if (mpfr_round_raw (MPFR_MANT (x), result, ysize_bits, + pstr->negative, precx, rnd, &res)) { /* overflow when rounding y */ MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT; @@ -769,12 +893,11 @@ exp ++; } - if (res == 0) /* fix ternary value */ - { - exact = exact && (pstr_size == pstr->prec); - if (!exact) - res = (pstr->negative) ? 1 : -1; - } + /* Note: if exact <> 0, then the approximation {result, ysize} is exact, + thus no double-rounding can occur: + (a) either the ternary value res is non-zero, and it is the correct + ternary value that we should return + (b) or the ternary value res is zero, and we should return 0. */ /* Set sign of x before exp since check_range needs a valid sign */ (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x); diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-12-31 11:18:05.621387341 +0000 +++ mpfr-4.0.1-b/src/version.c 2019-01-31 12:05:56.091705514 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p14"; + return "4.0.1-p15"; } diff -Naurd mpfr-4.0.1-a/tests/tstrtofr.c mpfr-4.0.1-b/tests/tstrtofr.c --- mpfr-4.0.1-a/tests/tstrtofr.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tstrtofr.c 2019-01-31 12:05:56.051705613 +0000 @@ -22,6 +22,12 @@ #include "mpfr-test.h" +/* The implicit \0 is useless, but we do not write num_to_text[62] otherwise + g++ complains. */ +static const char num_to_text36[] = "0123456789abcdefghijklmnopqrstuvwxyz"; +static const char num_to_text62[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" + "abcdefghijklmnopqrstuvwxyz"; + static void check_special (void) { @@ -898,6 +904,18 @@ { printf ("Check overflow failed (2) with:\n s='%s'\n x=", s); mpfr_dump (x); +#if defined(__GNUC__) + printf ("This failure is triggered by GCC bug 86554:\n" + " https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86554\n" + " https://gcc.gnu.org/bugzilla/show_bug.cgi?id=87276 " + "(about this test)\nWorkaround: disable code hoisting " + "with -fno-code-hoisting in CFLAGS.\n"); + /* Note: In Debian, this error is obtained with gcc-snapshot from + 20180908-1 to 20181127-1. With gcc-snapshot from 20181209-1 to + 20190102-1 (at least), the MPFR build no longer seems affected + in general, but using --with-gmp-build=... together with + --enable-assert still triggers this failure. */ +#endif exit (1); } mpfr_strtofr (x, "123456789E170141183460469231731687303715884105728", @@ -938,6 +956,20 @@ mpfr_dump (x); exit (1); } + mpfr_strtofr (x, "1@2305843009213693951", &s, 16, MPFR_RNDN); + if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x)) + { + printf ("Check overflow failed (8) with:\n s=%s\n x=", s); + mpfr_dump (x); + exit (1); + } + mpfr_strtofr (x, "1@2305843009213693951", &s, 17, MPFR_RNDN); + if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x)) + { + printf ("Check overflow failed (9) with:\n s=%s\n x=", s); + mpfr_dump (x); + exit (1); + } /* Check underflow */ mpfr_strtofr (x, "123456789E-2147483646", &s, 0, MPFR_RNDN); @@ -969,6 +1001,13 @@ mpfr_dump (x); exit (1); } + mpfr_strtofr (x, "1@-2305843009213693952", &s, 16, MPFR_RNDN); + if (s[0] != 0 || !MPFR_IS_ZERO (x) || !MPFR_IS_POS (x) ) + { + printf ("Check underflow failed (8) with:\n s='%s'\n x=", s); + mpfr_dump (x); + exit (1); + } mpfr_clear (x); } @@ -1238,12 +1277,225 @@ int inex; emin = mpfr_get_emin (); - mpfr_set_emin (-1073); - mpfr_set_emin (emin); mpfr_init2 (z, 53); + mpfr_set_emin (-1073); + /* with emin = -1073, the smallest positive number is 0.5*2^emin = 2^(-1074), + thus str should be rounded to 2^(-1074) with inex > 0 */ + inex = mpfr_strtofr (z, str, NULL, 10, MPFR_RNDN); + MPFR_ASSERTN(inex > 0 && mpfr_cmp_ui_2exp (z, 1, -1074) == 0); + mpfr_set_emin (-1074); + /* with emin = -1074, str should be rounded to 2^(-1075) with inex < 0 */ inex = mpfr_strtofr (z, str, NULL, 10, MPFR_RNDN); MPFR_ASSERTN(inex < 0 && mpfr_cmp_ui_2exp (z, 1, -1075) == 0); mpfr_clear (z); + mpfr_set_emin (emin); +} + +/* r13299 fails with 8-bit limbs (micro-gmp/8). */ +static void +bug20181127 (void) +{ + char s[] = "9.Z6nrLVSMG1bDcCF2ONJdX@-183295525"; /* base 58 */ + mpfr_t x, y; + + mpfr_inits2 (6, x, y, (mpfr_ptr) 0); + mpfr_set_ui_2exp (x, 5, -1073741701, MPFR_RNDN); + mpfr_strtofr (y, s, NULL, 58, MPFR_RNDZ); + if (! mpfr_equal_p (x, y)) + { + printf ("Error in bug20181127 on %s (base 58)\n", s); + printf ("Expected x = "); + mpfr_dump (x); + printf ("Got y = "); + mpfr_dump (y); + printf ("*Note* In base 58, x ~= "); + mpfr_out_str (stdout, 58, 8, x, MPFR_RNDN); + printf ("\n"); + exit (1); + } + mpfr_clears (x, y, (mpfr_ptr) 0); +} + +static void +coverage (void) +{ +#if _MPFR_EXP_FORMAT >= 3 && _MPFR_PREC_FORMAT == 3 && MPFR_PREC_BITS == 64 + char str3[] = "1@-2112009130072406892"; + char str31[] = "1@-511170973314085831"; + mpfr_t x; + int inex; + mpfr_exp_t emin; + + /* exercise assertion cy == 0 around line 698 of strtofr.c */ + emin = mpfr_get_emin (); + mpfr_set_emin (mpfr_get_emin_min ()); + /* emin = -4611686018427387903 on a 64-bit machine */ + mpfr_init2 (x, 1); + inex = mpfr_strtofr (x, str3, NULL, 3, MPFR_RNDN); + /* 3^-2112009130072406892 is slightly larger than (2^64)^-52303988630398057 + thus we should get inex < 0 */ + MPFR_ASSERTN(inex < 0); + MPFR_ASSERTN(mpfr_cmp_ui_2exp (x, 1, -52303988630398057 * 64) == 0); + inex = mpfr_strtofr (x, str31, NULL, 31, MPFR_RNDN); + /* 31^-511170973314085831 is slightly smaller than (2^64)^-39569396093273623 + thus we should get inex > 0 */ + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui_2exp (x, 1, -39569396093273623 * 64) == 0); + mpfr_clear (x); + mpfr_set_emin (emin); +#endif +} + +#define BSIZE 512 + +static void +random_tests (void) +{ + char s0[BSIZE], s1[BSIZE], s2[BSIZE+64]; + mpfr_t x0, x1, x2; + int prec, i; + + for (prec = MPFR_PREC_MIN; prec < 300; prec++) + { + mpfr_inits2 (prec, x0, x1, x2, (mpfr_ptr) 0); + + for (i = 0; i < 20; i++) + { + const char *num_to_text; + mpfr_exp_t e0, e1; + int base, j, k, neg; + int noteq = 0; + char d; + + /* We want the same exponent for x0 and its successor x1. + This is not possible for precision 1 in base 2. */ + do + base = 2 + (randlimb () % 61); + while (prec == 1 && base == 2); + + num_to_text = base <= 36 ? num_to_text36 : num_to_text62; + + do + { + /* Let's consider only positive numbers. We should test + negative numbers, but they should be built later, just + for the test itself. */ + tests_default_random (x0, 0, + mpfr_get_emin (), mpfr_get_emax (), 1); + mpfr_set (x1, x0, MPFR_RNDN); + mpfr_nextabove (x1); + mpfr_get_str (s0, &e0, base, BSIZE - 1, x0, MPFR_RNDU); + mpfr_get_str (s1, &e1, base, BSIZE - 1, x1, MPFR_RNDD); + } + while (! (mpfr_regular_p (x0) && mpfr_regular_p (x1) && e0 == e1)); + + /* 0 < x0 <= (s0,e) <= (s1,e) <= x1 with e = e0 = e1. + Let's build a string s2 randomly formed by: + - the common prefix of s0 and s1; + - some of the following digits of s0 (possibly none); + - the next digit of s0 + 1; + - some of the following digits of s1 (possibly none). + Then 0 < x0 <= (s0,e) < (s2,e) <= (s1,e) <= x1, and with + a very high probability that (s2,e) < (s1,e); noteq is + set to true in this case. + For instance, if: + s0 = 123456789 + s1 = 124012345 + one can have, e.g.: + s2 = 12345734 + s2 = 123556789 + s2 = 124 + s2 = 124012 + s2 is not taken completely randomly between s0 and s1, but it + will be built rather easily, and with the randomness of x0, + we should cover all cases, with s2 very close to s0, s2 very + close to s1, or not too close to either. */ + + neg = randlimb () & 1; + s2[0] = neg ? '-' : '+'; + s2[1] = '.'; + + for (j = 0; + MPFR_ASSERTN (s0[j] != 0 && s1[j] != 0), s0[j] == s1[j]; + j++) + s2[j+2] = s0[j]; + + /* k is the position of the first differing digit. */ + k = j; + + while (j < BSIZE - 2 && randlimb () % 8 != 0) + { + MPFR_ASSERTN (s0[j] != 0); + s2[j+2] = s0[j]; + j++; + } + + MPFR_ASSERTN (s0[j] != 0); + /* We will increment the next digit. Thus while s0[j] is the + maximum digit, go back until this is no longer the case + (the first digit after the common prefix cannot be the + maximum digit, so that we will stop early enough). */ + while ((d = s0[j]) == num_to_text[base - 1]) + j--; + noteq = j != k; + s2[j+2] = d = *(strchr (num_to_text, d) + 1); + if (d != s1[j]) + noteq = 1; + j++; + + while (j < BSIZE - 1 && randlimb () % 8 != 0) + { + MPFR_ASSERTN (s1[j] != 0); + s2[j+2] = s1[j]; + j++; + } + + sprintf (s2 + (j+2), "@%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e0); + + while (noteq == 0 && j < BSIZE - 1) + { + if (s1[j] != '0') + noteq = 1; + j++; + } + + if (neg) + { + mpfr_neg (x0, x0, MPFR_RNDN); + mpfr_neg (x1, x1, MPFR_RNDN); + } + + if (noteq) + { + mpfr_strtofr (x2, s2, NULL, base, MPFR_RNDZ); + if (! mpfr_equal_p (x2, x0)) + { + printf ("Error in random_tests for prec=%d i=%d base=%d\n", + prec, i, base); + printf ("s0 = %s\ns1 = %s\ns2 = %s\n", s0, s1, s2); + printf ("x0 = "); + mpfr_dump (x0); + printf ("x2 = "); + mpfr_dump (x2); + exit (1); + } + } + + mpfr_strtofr (x2, s2, NULL, base, MPFR_RNDA); + if (! mpfr_equal_p (x2, x1)) + { + printf ("Error in random_tests for prec=%d i=%d base=%d\n", + prec, i, base); + printf ("s0 = %s\ns1 = %s\ns2 = %s\n", s0, s1, s2); + printf ("x1 = "); + mpfr_dump (x1); + printf ("x2 = "); + mpfr_dump (x2); + exit (1); + } + } + mpfr_clears (x0, x1, x2, (mpfr_ptr) 0); + } } int @@ -1251,6 +1503,7 @@ { tests_start_mpfr (); + coverage (); check_special(); check_reftable (); check_parse (); @@ -1262,6 +1515,8 @@ bug20120829 (); bug20161217 (); bug20170308 (); + bug20181127 (); + random_tests (); tests_end_mpfr (); return 0;