[src/strtofr.c] Fixed various issues. In particular, the error analysis

  with associated code was incorrect (due to the bad correction r8384
  of a past bug, later really fixed in r11056). Also adapted the code
  to work with small-size limbs.
[tests/tstrtofr.c] Added tests, including random tests. In particular,
  some part of the code was tested only on hard-to-round cases, meaning
  that some potential issues could not be detected (intermediate results
  were thrown out due to the next Ziv iteration). Moreover, in case of
  failure of some particular test, output the probable cause of this
  failure (GCC bug 86554) with a workaround.
(merged changesets associated with these files r12566-13306,13364-13365
from the trunk)
Note: There may still remain issues with mpfr_strtofr, to be checked
later. At least, the tests do not fail.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/4.0@13366 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 555 insertions, 177 deletions

diff --git a/src/strtofr.c b/src/strtofr.c
index 7dc01eba6..f2a8964b4 100644
--- a/src/strtofr.c
+++ b/src/strtofr.c
@@ -36,79 +36,80 @@ struct parsed_string {
                               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 @@ parse_string (mpfr_t x, struct parsed_string *pstr,
   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 @@ parse_string (mpfr_t x, struct parsed_string *pstr,
 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));
-
-  /* compute the value y of the leading characters as long as rounding is not
-     possible */
+  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.
+     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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
           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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
           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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
 
           /* 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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
           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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
                               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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
                               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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
       /* 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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
 
       /* 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 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
       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 --git a/tests/tstrtofr.c b/tests/tstrtofr.c
index a6efac63d..d53b231a2 100644
--- a/tests/tstrtofr.c
+++ b/tests/tstrtofr.c
@@ -22,6 +22,12 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 
 #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 @@ check_overflow (void)
     {
       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 @@ check_overflow (void)
       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 @@ check_overflow (void)
       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 @@ bug20170308 (void)
   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 @@ main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
+  coverage ();
   check_special();
   check_reftable ();
   check_parse ();
@@ -1262,6 +1515,8 @@ main (int argc, char *argv[])
   bug20120829 ();
   bug20161217 ();
   bug20170308 ();
+  bug20181127 ();
+  random_tests ();
 
   tests_end_mpfr ();
   return 0;