The fix for bug 6604 "incorrect directed rounding in mpfr_strtofr"2.3

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

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