Fixed various bugs in mpfr_can_round_raw, which affected mpfr_can_round:

there could still be false positives (i.e. mpfr_can_round could say that
rounding was possible while correct rounding was not guaranteed), and
also false negatives, some of which could yield infinite Ziv loops in
user code in practice.
Added tests triggering these bugs, in particular a comprehensive test
against a naive implementation.
(merged changesets
  r10679-10686,10717-10718,10743,10746-10748,10752,10754,10756
from the trunk)

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

5 files changed, 393 insertions, 62 deletions

diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h
index 1999b5abc..e9c0db62f 100644
--- a/src/mpfr-impl.h
+++ b/src/mpfr-impl.h
@@ -1889,6 +1889,7 @@ __MPFR_DECLSPEC mpfr_exp_t mpfr_ceil_mul _MPFR_PROTO ((mpfr_exp_t, int, int));
 __MPFR_DECLSPEC int mpfr_exp_2 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t));
 __MPFR_DECLSPEC int mpfr_exp_3 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t));
 __MPFR_DECLSPEC int mpfr_powerof2_raw _MPFR_PROTO ((mpfr_srcptr));
+__MPFR_DECLSPEC int mpfr_powerof2_raw2 (const mp_limb_t *, mp_size_t);
 
 __MPFR_DECLSPEC int mpfr_pow_general _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
                            mpfr_srcptr, mpfr_rnd_t, int, mpfr_save_expo_t *));
diff --git a/src/powerof2.c b/src/powerof2.c
index f10ace0ad..53f0551cd 100644
--- a/src/powerof2.c
+++ b/src/powerof2.c
@@ -31,16 +31,17 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 int
 mpfr_powerof2_raw (mpfr_srcptr x)
 {
-  mp_limb_t *xp;
-  mp_size_t xn;
-
   /* This is an internal function, and we may call it with some
      wrong numbers (ie good mantissa but wrong flags or exp)
      So we don't want to test if it is a pure FP.
      MPFR_ASSERTN(MPFR_IS_PURE_FP(x)); */
-  xp = MPFR_MANT(x);
-  xn = (MPFR_PREC(x) - 1) / GMP_NUMB_BITS;
-  if (xp[xn] != MPFR_LIMB_HIGHBIT)
+  return mpfr_powerof2_raw2 (MPFR_MANT(x), MPFR_LIMB_SIZE(x));
+}
+
+int
+mpfr_powerof2_raw2 (const mp_limb_t *xp, mp_size_t xn)
+{
+  if (xp[--xn] != MPFR_LIMB_HIGHBIT)
     return 0;
   while (xn > 0)
     if (xp[--xn] != 0)
diff --git a/src/round_p.c b/src/round_p.c
index f19a66433..8e2c0bebd 100644
--- a/src/round_p.c
+++ b/src/round_p.c
@@ -35,10 +35,16 @@ mpfr_round_p (mp_limb_t *bp, mp_size_t bn, mpfr_exp_t err0, mpfr_prec_t prec)
 
   i1 = mpfr_round_p_2 (bp, bn, err0, prec);
 
-  /* compare with mpfr_can_round_raw */
+  /* Note: since revision 10747, mpfr_can_round_raw is supposed to be always
+     correct, whereas mpfr_round_p_2 might return 0 in some cases where one
+     could round, for example with err0=67 and prec=54:
+     b = 1111101101010001100011111011100010100011101111011011101111111111
+     thus we cannot compare i1 and i2, we only can check that we don't have
+     i1 <> 0 and i2 = 0.
+  */
   i2 = mpfr_can_round_raw (bp, bn, MPFR_SIGN_POS, err0,
                            MPFR_RNDN, MPFR_RNDZ, prec);
-  if (i1 != i2)
+  if (i1 && (i2 == 0))
     {
       fprintf (stderr, "mpfr_round_p(%d) != mpfr_can_round(%d)!\n"
                "bn = %lu, err0 = %ld, prec = %lu\nbp = ", i1, i2,
diff --git a/src/round_prec.c b/src/round_prec.c
index c99d63e55..e5512b439 100644
--- a/src/round_prec.c
+++ b/src/round_prec.c
@@ -138,47 +138,173 @@ mpfr_can_round (mpfr_srcptr b, mpfr_exp_t err, mpfr_rnd_t rnd1,
 }
 
 int
-mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
+mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
                     mpfr_rnd_t rnd1, mpfr_rnd_t rnd2, mpfr_prec_t prec)
 {
-  mpfr_prec_t err, prec0 = prec;
+  mpfr_prec_t prec2;
   mp_size_t k, k1, tn;
   int s, s1;
   mp_limb_t cc, cc2;
   mp_limb_t *tmp;
+  mp_limb_t cy = 0, tmp_hi;
+  int res;
   MPFR_TMP_DECL(marker);
 
-  MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT);
+  /* Since mpfr_can_round is a function in the API, use MPFR_ASSERTN.
+     The specification makes sense only for prec >= 1. */
+  MPFR_ASSERTN (prec >= 1);
 
-  if (MPFR_UNLIKELY(err0 < 0 || (mpfr_uexp_t) err0 <= prec))
-    return 0;  /* can't round */
+  MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT);
 
   MPFR_ASSERT_SIGN(neg);
   neg = MPFR_IS_NEG_SIGN(neg);
 
   /* Transform RNDD and RNDU to Zero / Away */
-  MPFR_ASSERTD((neg == 0) || (neg == 1));
+  MPFR_ASSERTD (neg == 0 || neg == 1);
   if (rnd1 != MPFR_RNDN)
     rnd1 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDZ : MPFR_RNDA;
   if (rnd2 != MPFR_RNDN)
     rnd2 = MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDZ : MPFR_RNDA;
 
+  /* For err < prec (+1 for rnd1=RNDN), we can never round correctly, since
+     the error is at least 2*ulp(b) >= ulp(round(b)).
+     However for err = prec (+1 for rnd1=RNDN), we can round correctly in some
+     rare cases where ulp(b) = 1/2*ulp(U) [see below for the definition of U],
+     which implies rnd1 = RNDZ or RNDN, and rnd2 = RNDA or RNDN. */
+
+  if (MPFR_UNLIKELY (err < prec + (rnd1 == MPFR_RNDN) ||
+                     (err == prec + (rnd1 == MPFR_RNDN) &&
+                      (rnd1 == MPFR_RNDA ||
+                       rnd2 == MPFR_RNDZ))))
+    return 0;  /* can't round */
+
+  /* As a consequence... */
+  MPFR_ASSERTD (err >= prec);
+
+  /* The bound c on the error |x-b| is: c = 2^(MPFR_EXP(b)-err) <= b/2.
+   * So, we now know that x and b have the same sign. By symmetry,
+   * assume x > 0 and b > 0. We have: L <= x <= U, where, depending
+   * on rnd1:
+   *   MPFR_RNDN: L = b-c, U = b+c
+   *   MPFR_RNDZ: L = b,   U = b+c
+   *   MPFR_RNDA: L = b-c, U = b
+   *
+   * We can round x iff round(L,prec,rnd2) = round(U,prec,rnd2).
+   */
+
   if (MPFR_UNLIKELY (prec > (mpfr_prec_t) bn * GMP_NUMB_BITS))
-    { /* Then prec < PREC(b): we can round:
-         (i) in rounding to the nearest iff err0 >= prec + 2
+    { /* Then prec > PREC(b): we can round:
+         (i) in rounding to the nearest as long as err >= prec + 2.
+             When err = prec + 1 and b is not a power
+             of two (so that a change of binade cannot occur), then one
+             can round to nearest thanks to the even rounding rule (in the
+             target precision prec, the significand of b ends with a 0).
+             When err = prec + 1 and b is a power of two, when rnd1 = RNDZ one
+             can round too.
          (ii) in directed rounding mode iff rnd1 is compatible with rnd2
-              and err0 >= prec + 1, unless b = 2^k and rnd1=rnd2=RNDA in
-              which case we need err0 >= prec + 2. */
+              and err >= prec + 1, unless b = 2^k and rnd1 = RNDA or RNDN in
+              which case we need err >= prec + 2.
+      */
+      if ((rnd1 == rnd2 || rnd2 == MPFR_RNDN) && err >= prec + 1)
+        {
+          if (rnd1 != MPFR_RNDZ &&
+              err == prec + 1 &&
+              mpfr_powerof2_raw2 (bp, bn))
+            return 0;
+          else
+            return 1;
+        }
+      return 0;
+    }
+
+  /* now prec <= bn * GMP_NUMB_BITS */
+
+  if (MPFR_UNLIKELY (err > (mpfr_prec_t) bn * GMP_NUMB_BITS))
+    {
+      /* we distinguish the case where b is a power of two:
+         rnd1 rnd2 can round?
+         RNDZ RNDZ ok
+         RNDZ RNDA no
+         RNDZ RNDN ok
+         RNDA RNDZ no
+         RNDA RNDA ok except when err = prec + 1
+         RNDA RNDN ok except when err = prec + 1
+         RNDN RNDZ no
+         RNDN RNDA no
+         RNDN RNDN ok except when err = prec + 1 */
+      if (mpfr_powerof2_raw2 (bp, bn))
+        {
+          if ((rnd2 == MPFR_RNDZ || rnd2 == MPFR_RNDA) && rnd1 != rnd2)
+            return 0;
+          else if (rnd1 == MPFR_RNDZ)
+            return 1; /* RNDZ RNDZ and RNDZ RNDN */
+          else
+            return err > prec + 1;
+        }
+
+      /* now the general case where b is not a power of two:
+         rnd1 rnd2 can round?
+         RNDZ RNDZ ok
+         RNDZ RNDA except when b is representable in precision 'prec'
+         RNDZ RNDN except when b is the middle of two representable numbers in
+                   precision 'prec' and b ends with 'xxx0[1]',
+                   or b is representable in precision 'prec'
+                   and err = prec + 1 and b ends with '1'.
+         RNDA RNDZ except when b is representable in precision 'prec'
+         RNDA RNDA ok
+         RNDA RNDN except when b is the middle of two representable numbers in
+                   precision 'prec' and b ends with 'xxx1[1]',
+                   or b is representable in precision 'prec'
+                   and err = prec + 1 and b ends with '1'.
+         RNDN RNDZ except when b is representable in precision 'prec'
+         RNDN RNDA except when b is representable in precision 'prec'
+         RNDN RNDN except when b is the middle of two representable numbers in
+                   precision 'prec', or b is representable in precision 'prec'
+                   and err = prec + 1 and b ends with '1'. */
       if (rnd2 == MPFR_RNDN)
-        return (mpfr_uexp_t) err0 - 2 >= prec;
+        {
+          if (err == prec + 1 && (bp[0] & 1))
+            return 0; /* err == prec + 1 implies prec = bn * GMP_NUMB_BITS */
+          if (prec < (mpfr_prec_t) bn * GMP_NUMB_BITS)
+            {
+              k1 = MPFR_PREC2LIMBS (prec + 1);
+              MPFR_UNSIGNED_MINUS_MODULO(s1, prec + 1);
+              if (((bp[bn - k1] >> s1) & 1) &&
+                  mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec + 1) == 0)
+                { /* b is the middle of two representable numbers */
+                  if (rnd1 == MPFR_RNDN)
+                    return 0;
+                  k1 = MPFR_PREC2LIMBS (prec);
+                  MPFR_UNSIGNED_MINUS_MODULO(s1, prec);
+                  return (rnd1 == MPFR_RNDZ) ^
+                    (((bp[bn - k1] >> s1) & 1) == 0);
+                }
+            }
+          return 1;
+        }
+      else if (rnd1 == rnd2)
+        {
+          if (rnd1 == MPFR_RNDN && prec < (mpfr_prec_t) bn * GMP_NUMB_BITS)
+            {
+              /* then rnd2 = RNDN, and for prec = bn * GMP_NUMB_BITS we cannot
+                 have b the middle of two representable numbers */
+              k1 = MPFR_PREC2LIMBS (prec + 1);
+              MPFR_UNSIGNED_MINUS_MODULO(s1, prec + 1);
+              if (((bp[bn - k1] >> s1) & 1) &&
+                  mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec + 1) == 0)
+                /* b is representable in precision prec+1 and ends with a 1 */
+                return 0;
+              else
+                return 1;
+            }
+          else
+            return 1;
+        }
       else
-        return (rnd1 == rnd2) && (mpfr_uexp_t) err0 - 2 >= prec;
+        return mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec) != 0;
     }
 
-  /* if the error is smaller than ulp(b), then anyway it will propagate
-     up to ulp(b) */
-  err = ((mpfr_uexp_t) err0 > (mpfr_prec_t) bn * GMP_NUMB_BITS) ?
-    (mpfr_prec_t) bn * GMP_NUMB_BITS : (mpfr_prec_t) err0;
+  /* now err <= bn * GMP_NUMB_BITS */
 
   /* warning: if k = m*GMP_NUMB_BITS, consider limb m-1 and not m */
   k = (err - 1) / GMP_NUMB_BITS;
@@ -196,7 +322,7 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
      Warning! The number with updated bn may no longer be normalized. */
   k -= k1;
   bn -= k1;
-  prec -= (mpfr_prec_t) k1 * GMP_NUMB_BITS;
+  prec2 = prec - (mpfr_prec_t) k1 * GMP_NUMB_BITS;
 
   /* We can decide of the correct rounding if rnd2(b-eps) and rnd2(b+eps)
      give the same result to the target precision 'prec', i.e., if when
@@ -215,64 +341,150 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
   switch (rnd1)
     {
     case MPFR_RNDZ:
-      /* Round to Zero */
+      /* rnd1 = Round to Zero */
       cc = (bp[bn - 1] >> s1) & 1;
       /* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec),
          and 0 otherwise */
-      cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
+      cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2);
       /* cc is the new value of bit s1 in bp[bn-1] after rounding 'rnd2' */
 
       /* now round b + 2^(MPFR_EXP(b)-err) */
-      mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
-      /* if there was a carry here, then necessarily bit s1 of bp[bn-1]
-         changed, thus we surely cannot round for directed rounding, but this
-         will be detected below, with cc2 != cc */
+      cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+      /* propagate carry up to most significant limb */
+      for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
+        cy = ~bp[bn + tn] == 0;
+      if (cy == 0 && err == prec)
+        {
+          res = 0;
+          goto end;
+        }
+      if (MPFR_UNLIKELY(cy))
+        {
+          /* when a carry occurs, we have b < 2^h <= b+c, we can round iff:
+             rnd2 = RNDZ: never, since b and b+c round to different values;
+             rnd2 = RNDA: when b+c is an exact power of two, and err > prec
+                          (since for err = prec, b = 2^h - 1/2*ulp(2^h) is
+                          exactly representable and thus rounds to itself);
+             rnd2 = RNDN: whenever cc = 0, since err >= prec implies
+                          c <= ulp(b) = 1/2*ulp(2^h), thus b+c rounds to 2^h,
+                          and b+c >= 2^h implies that bit 'prec' of b is 1,
+                          thus cc = 0 means that b is rounded to 2^h too. */
+          res = (rnd2 == MPFR_RNDZ) ? 0
+            : (rnd2 == MPFR_RNDA) ? (err > prec && k == bn && tmp[0] == 0)
+            : cc == 0;
+          goto end;
+        }
       break;
     case MPFR_RNDN:
-      /* Round to nearest */
+      /* rnd1 = Round to nearest */
 
       /* first round b+2^(MPFR_EXP(b)-err) */
-      mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
-      /* same remark as above in case a carry occurs in mpn_add_1() */
+      cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+      /* propagate carry up to most significant limb */
+      for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
+        cy = ~bp[bn + tn] == 0;
       cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
-      cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
+      cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2);
       /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */
-
+      if (MPFR_UNLIKELY (cy != 0))
+        {
+          /* when a carry occurs, we have b-c < b < 2^h <= b+c, we can round
+             iff:
+             rnd2 = RNDZ: never, since b-c and b+c round to different values;
+             rnd2 = RNDA: when b+c is an exact power of two, and
+                          err > prec + 1 (since for err <= prec + 1,
+                          b-c <= 2^h - 1/2*ulp(2^h) is exactly representable
+                          and thus rounds to itself);
+             rnd2 = RNDN: whenever err > prec + 1, since for err = prec + 1,
+                          b+c rounds to 2^h, and b-c rounds to nextbelow(2^h).
+                          For err > prec + 1, c <= 1/4*ulp(b) <= 1/8*ulp(2^h),
+                          thus
+                          2^h - 1/4*ulp(b) <= b-c < b+c <= 2^h + 1/8*ulp(2^h),
+                          therefore both b-c and b+c round to 2^h. */
+          res = (rnd2 == MPFR_RNDZ) ? 0
+            : (rnd2 == MPFR_RNDA) ? (err > prec + 1 && k == bn && tmp[0] == 0)
+            : err > prec + 1;
+          goto end;
+        }
     subtract_eps:
-      /* now round b-2^(MPFR_EXP(b)-err) */
-      cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+      /* now round b-2^(MPFR_EXP(b)-err), this happens for
+         rnd1 = RNDN or RNDA */
+      MPFR_ASSERTD(rnd1 == MPFR_RNDN || rnd1 == MPFR_RNDA);
+      cy = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
       /* propagate the potential borrow up to the most significant limb
          (it cannot propagate further since the most significant limb is
-         at least MPFR_LIMB_HIGHBIT) */
-      for (tn = 0; tn + 1 < k1 && (cc2 != 0); tn ++)
-        cc2 = bp[bn + tn] == 0;
-      /* We have an exponent decrease when either:
-           (i) k1 = 0 and tmp[bn-1] < MPFR_LIMB_HIGHBIT
-           (ii) k1 > 0 and cc <> 0 and bp[bn + tn] = MPFR_LIMB_HIGHBIT
-                (then necessarily tn = k1-1).
-         Then for directed rounding we cannot round,
-         and for rounding to nearest we cannot round when err = prec + 1.
+         at least MPFR_LIMB_HIGHBIT).
+         Note: we use the same limb tmp[bn-1] to subtract. */
+      tmp_hi = tmp[bn - 1];
+      for (tn = 0; tn < k1 && cy != 0; tn ++)
+        cy = mpn_sub_1 (&tmp_hi, bp + bn + tn, 1, cy);
+      /* We have an exponent decrease when tn = k1 and
+         tmp[bn-1] < MPFR_LIMB_HIGHBIT:
+         b-c < 2^h <= b (for RNDA) or b+c (for RNDN).
+         Then we surely cannot round when rnd2 = RNDZ, since b or b+c round to
+         a value >= 2^h, and b-c rounds to a value < 2^h.
+         We also surely cannot round when (rnd1,rnd2) = (RNDN,RNDA), since
+         b-c rounds to a value <= 2^h, and b+c > 2^h rounds to a value > 2^h.
+         It thus remains:
+         (rnd1,rnd2) = (RNDA,RNDA), (RNDA,RNDN) and (RNDN,RNDN).
+         For (RNDA,RNDA) we can round only when b-c and b round to 2^h, which
+         implies b = 2^h and err > prec (which is true in that case):
+         a necessary condition is that cc = 0.
+         For (RNDA,RNDN) we can round only when b-c and b round to 2^h, which
+         implies b-c >= 2^h - 1/4*ulp(2^h), and b <= 2^h + 1/2*ulp(2^h);
+         since ulp(2^h) = ulp(b), this implies c <= 3/4*ulp(b), thus
+         err > prec.
+         For (RNDN,RNDN) we can round only when b-c and b+c round to 2^h,
+         which implies b-c >= 2^h - 1/4*ulp(2^h), and
+         b+c <= 2^h + 1/2*ulp(2^h);
+         since ulp(2^h) = ulp(b), this implies 2*c <= 3/4*ulp(b), thus
+         err > prec+1.
       */
-      if (((k1 == 0 && tmp[bn - 1] < MPFR_LIMB_HIGHBIT) ||
-           (k1 != 0 && cc2 != 0 && bp[bn + tn] == MPFR_LIMB_HIGHBIT)) &&
-          (rnd2 != MPFR_RNDN || err0 == prec0 + 1))
+      if (tn == k1 && tmp_hi < MPFR_LIMB_HIGHBIT) /* exponent decrease */
+        {
+          if (rnd2 == MPFR_RNDZ || (rnd1 == MPFR_RNDN && rnd2 == MPFR_RNDA) ||
+              cc != 0 /* b or b+c does not round to 2^h */)
+            {
+              res = 0;
+              goto end;
+            }
+          /* in that case since the most significant bit of tmp is 0, we
+             should consider one more bit; res = 0 when b-c does not round
+             to 2^h. */
+          res = mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2 + 1) != 0;
+          goto end;
+        }
+      if (err == prec + (rnd1 == MPFR_RNDN))
         {
-          MPFR_TMP_FREE(marker);
-          return 0;
+          /* No exponent increase nor decrease, thus we have |U-L| = ulp(b).
+             For rnd2 = RNDZ or RNDA, either [L,U] contains one representable
+             number in the target precision, and then L and U round
+             differently; or both L and U are representable: they round
+             differently too; thus in all cases we cannot round.
+             For rnd2 = RNDN, the only case where we can round is when the
+             middle of [L,U] (i.e. b) is representable, and ends with a 0. */
+          res = (rnd2 == MPFR_RNDN && (((bp[bn - 1] >> s1) & 1) == 0) &&
+                 mpfr_round_raw2 (bp, bn, neg, MPFR_RNDZ, prec2) ==
+                 mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec2));
+          goto end;
         }
       break;
     default:
-      /* Round away */
+      /* rnd1 = Round away */
+      MPFR_ASSERTD (rnd1 == MPFR_RNDA);
       cc = (bp[bn - 1] >> s1) & 1;
-      cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
+      /* the mpfr_round_raw2() call below returns whether one should add 1 or
+         not for rounding */
+      cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2);
       /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */
 
       goto subtract_eps;
     }
 
   cc2 = (tmp[bn - 1] >> s1) & 1;
-  cc2 ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
+  res = cc == (cc2 ^ mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2));
 
+ end:
   MPFR_TMP_FREE(marker);
-  return cc == cc2;
+  return res;
 }
diff --git a/tests/tcan_round.c b/tests/tcan_round.c
index a415d2f5d..108af8862 100644
--- a/tests/tcan_round.c
+++ b/tests/tcan_round.c
@@ -25,6 +25,44 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 
 #include "mpfr-test.h"
 
+/* Simple cases where err - prec is large enough.
+   One can get failures as a consequence of r9883. */
+static void
+test_simple (void)
+{
+  int t[4] = { 2, 3, -2, -3 };  /* test powers of 2 and non-powers of 2 */
+  int i;
+  int r1, r2;
+
+  for (i = 0; i < 4; i++)
+    RND_LOOP (r1)
+      RND_LOOP (r2)
+        {
+          mpfr_t b;
+          int p, err, prec, inex, c;
+
+          p = 12 + (randlimb() % (2 * GMP_NUMB_BITS));
+          err = p - 3;
+          prec = 4;
+          mpfr_init2 (b, p);
+          inex = mpfr_set_si (b, t[i], MPFR_RNDN);
+          MPFR_ASSERTN (inex == 0);
+          c = mpfr_can_round (b, err, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, prec);
+          /* If r1 == r2, we can round.
+             TODO: complete this test for r1 != r2. */
+          if (r1 == r2 && !c)
+            {
+              printf ("Error in test_simple for i=%d,"
+                      " err=%d r1=%s, r2=%s, p=%d\n", i, err,
+                      mpfr_print_rnd_mode ((mpfr_rnd_t) r1),
+                      mpfr_print_rnd_mode ((mpfr_rnd_t) r2), p);
+              printf ("b="); mpfr_dump (b);
+              exit (1);
+            }
+          mpfr_clear (b);
+        }
+}
+
 #define MAX_LIMB_SIZE 100
 
 static void
@@ -87,12 +125,7 @@ test_pow2 (mpfr_exp_t i, mpfr_prec_t px, mpfr_rnd_t r1, mpfr_rnd_t r2,
     MPFR_IS_LIKE_RNDU (r1, MPFR_SIGN_POS) ?
     (MPFR_IS_LIKE_RNDD (r2, MPFR_SIGN_POS) ? 0 : prec < i) :
     (r2 != MPFR_RNDN ? 0 : prec < i);
-  /* We only require mpfr_can_round to return 1 when we can really
-     round, it is allowed to return 0 in some rare boundary cases,
-     for example when x = 2^k and the error is 0.25 ulp.
-     Note: if this changes in the future, the test could be improved by
-     removing the "&& expected_b == 0" below. */
-  if (b != expected_b && expected_b == 0)
+  if (b != expected_b)
     {
       printf ("Error for x=2^%d, px=%lu, err=%d, r1=%s, r2=%s, prec=%d\n",
               (int) i, (unsigned long) px, (int) i + 1,
@@ -118,6 +151,80 @@ test_pow2 (mpfr_exp_t i, mpfr_prec_t px, mpfr_rnd_t r1, mpfr_rnd_t r2,
   mpfr_clear (x);
 }
 
+static void
+check_can_round (void)
+{
+  mpfr_t x, xinf, xsup, yinf, ysup;
+  int precx, precy, err;
+  int rnd1, rnd2;
+  int i, u[3] = { 0, 1, 256 };
+  int inex;
+  int expected, got;
+
+  mpfr_inits2 (4 * GMP_NUMB_BITS, x, xinf, xsup, yinf, ysup, (mpfr_ptr) 0);
+
+  for (precx = 3 * GMP_NUMB_BITS - 3; precx <= 3 * GMP_NUMB_BITS + 3; precx++)
+    {
+      mpfr_set_prec (x, precx);
+      for (precy = precx - 4; precy <= precx + 4; precy++)
+        {
+          mpfr_set_prec (yinf, precy);
+          mpfr_set_prec (ysup, precy);
+
+          for (i = 0; i <= 3; i++)
+            {
+              if (i <= 2)
+                {
+                  /* Test x = 2^(precx - 1) + u */
+                  mpfr_set_ui_2exp (x, 1, precx - 1, MPFR_RNDN);
+                  mpfr_add_ui (x, x, u[i], MPFR_RNDN);
+                }
+              else
+                {
+                  /* Test x = 2^precx - 1 */
+                  mpfr_set_ui_2exp (x, 1, precx, MPFR_RNDN);
+                  mpfr_sub_ui (x, x, 1, MPFR_RNDN);
+                }
+              MPFR_ASSERTN (mpfr_get_exp (x) == precx);
+
+              for (err = precy; err <= precy + 3; err++)
+                {
+                  mpfr_set_ui_2exp (xinf, 1, precx - err, MPFR_RNDN);
+                  inex = mpfr_add (xsup, x, xinf, MPFR_RNDN);
+                  MPFR_ASSERTN (inex == 0);
+                  inex = mpfr_sub (xinf, x, xinf, MPFR_RNDN);
+                  MPFR_ASSERTN (inex == 0);
+                  RND_LOOP (rnd1)
+                    RND_LOOP (rnd2)
+                      {
+                        mpfr_set (yinf, MPFR_IS_LIKE_RNDD (rnd1, 1) ?
+                                  x : xinf, (mpfr_rnd_t) rnd2);
+                        mpfr_set (ysup, MPFR_IS_LIKE_RNDU (rnd1, 1) ?
+                                  x : xsup, (mpfr_rnd_t) rnd2);
+                        expected = !! mpfr_equal_p (yinf, ysup);
+                        got = !! mpfr_can_round (x, err, (mpfr_rnd_t) rnd1,
+                                                 (mpfr_rnd_t) rnd2, precy);
+                        if (got != expected)
+                          {
+                            printf ("Error in check_can_round on:\n"
+                                    "precx=%d, precy=%d, i=%d, err=%d, "
+                                    "rnd1=%s, rnd2=%s: got %d\n",
+                                    precx, precy, i, err,
+                                    mpfr_print_rnd_mode ((mpfr_rnd_t) rnd1),
+                                    mpfr_print_rnd_mode ((mpfr_rnd_t) rnd2),
+                                    got);
+                            printf ("x="); mpfr_dump (x);
+                            exit (1);
+                          }
+                      }
+                }
+            }
+        }
+    }
+
+  mpfr_clears (x, xinf, xsup, yinf, ysup, (mpfr_ptr) 0);
+}
+
 int
 main (void)
 {
@@ -128,6 +235,8 @@ main (void)
 
   tests_start_mpfr ();
 
+  test_simple ();
+
   /* checks that rounds to nearest sets the last
      bit to zero in case of equal distance */
   mpfr_init2 (x, 59);
@@ -201,6 +310,8 @@ main (void)
 
   mpfr_clear (x);
 
+  check_can_round ();
+
   check_round_p ();
 
   tests_end_mpfr ();