Fixed bug in mpfr_can_round_raw, which affected mpfr_can_round: the

result could be true instead of false in case of a change of binade
(exponent decrease) on the approximation interval. At the same time,
make sure that the number is normalized, and ditto for mpfr_round_p;
otherwise the semantic is not clear. Thus mpfr_div, which could call
mpfr_round_p with an unnormalized number, had to be fixed.
(merged changesets r9881,9883-9890,9896-9904,9932,10027 from the trunk)

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

1 files changed, 89 insertions, 8 deletions

diff --git a/tests/tcan_round.c b/tests/tcan_round.c
index f114d282d..a415d2f5d 100644
--- a/tests/tcan_round.c
+++ b/tests/tcan_round.c
@@ -1,4 +1,4 @@
-/* Test file for mpfr_can_round.
+/* Test file for mpfr_can_round and mpfr_round_p.
 
 Copyright 1999, 2001-2016 Free Software Foundation, Inc.
 Contributed by the AriC and Caramba projects, INRIA.
@@ -41,6 +41,8 @@ check_round_p (void)
       /* avoid mpn_random which leaks memory */
       for (i = 0; i < n; i++)
         buf[i] = randlimb ();
+      /* force the number to be normalized */
+      buf[n - 1] |= MPFR_LIMB_HIGHBIT;
       p = randlimb() % ((n-1) * GMP_NUMB_BITS) + MPFR_PREC_MIN;
       err = p + randlimb () % GMP_NUMB_BITS;
       r1 = mpfr_round_p (buf, n, err, p);
@@ -57,11 +59,72 @@ check_round_p (void)
     }
 }
 
+/* check x=2^i with precision px, error at most 1, and target precision prec */
+static void
+test_pow2 (mpfr_exp_t i, mpfr_prec_t px, mpfr_rnd_t r1, mpfr_rnd_t r2,
+           mpfr_prec_t prec)
+{
+  mpfr_t x;
+  int b, expected_b, b2;
+
+  mpfr_init2 (x, px);
+  mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN);
+  b = !!mpfr_can_round (x, i+1, r1, r2, prec);
+  /* Note: If mpfr_can_round succeeds for both
+     (r1,r2) = (MPFR_RNDD,MPFR_RNDN) and
+     (r1,r2) = (MPFR_RNDU,MPFR_RNDN), then it should succeed for
+     (r1,r2) = (MPFR_RNDN,MPFR_RNDN). So, the condition on prec below
+     for r1 = MPFR_RNDN should be the most restrictive between those
+     for r1 = any directed rounding mode.
+     For r1 like MPFR_RNDA, the unrounded, unknown number may be anyone
+     in [2^i-1,i]. As both 2^i-1 and 2^i fit on i bits, one cannot round
+     in any precision >= i bits, hence the condition prec < i; prec = i-1
+     will work here for r2 = MPFR_RNDN thanks to the even-rounding rule
+     (and also with rounding ties away from zero). */
+  expected_b =
+    MPFR_IS_LIKE_RNDD (r1, MPFR_SIGN_POS) ?
+    (MPFR_IS_LIKE_RNDU (r2, MPFR_SIGN_POS) ? 0 : prec <= i) :
+    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)
+    {
+      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,
+              mpfr_print_rnd_mode (r1), mpfr_print_rnd_mode (r2), (int) prec);
+      printf ("Expected %d, got %d\n", expected_b, b);
+      exit (1);
+    }
+
+  if (r1 == MPFR_RNDN && r2 == MPFR_RNDZ)
+    {
+      /* Similar test to the one done in src/round_p.c
+         for MPFR_WANT_ASSERT >= 2. */
+      b2 = !!mpfr_round_p (MPFR_MANT(x), MPFR_LIMB_SIZE(x), i+1, prec);
+      if (b2 != b)
+        {
+          printf ("Error for x=2^%d, px=%lu, err=%d, prec=%d\n",
+                  (int) i, (unsigned long) px, (int) i + 1, (int) prec);
+          printf ("mpfr_can_round gave %d, mpfr_round_p gave %d\n", b, b2);
+          exit (1);
+        }
+    }
+
+  mpfr_clear (x);
+}
+
 int
 main (void)
 {
   mpfr_t x;
-  mpfr_prec_t i, j;
+  mpfr_prec_t i, j, k;
+  int r1, r2;
+  int n;
 
   tests_start_mpfr ();
 
@@ -111,12 +174,30 @@ main (void)
   mpfr_set_str (x, "0.ff4ca619c76ba69", 16, MPFR_RNDZ);
   for (i = 30; i < 99; i++)
     for (j = 30; j < 99; j++)
-      {
-        int r1, r2;
-        for (r1 = 0; r1 < MPFR_RND_MAX ; r1++)
-          for (r2 = 0; r2 < MPFR_RND_MAX ; r2++)
-            mpfr_can_round (x, i, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j); /* test for assertions */
-      }
+      for (r1 = 0; r1 < MPFR_RND_MAX; r1++)
+        for (r2 = 0; r2 < MPFR_RND_MAX; r2++)
+          {
+            /* test for assertions */
+            mpfr_can_round (x, i, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j);
+          }
+
+  test_pow2 (32, 32, MPFR_RNDN, MPFR_RNDN, 32);
+  test_pow2 (174, 174, MPFR_RNDN, MPFR_RNDN, 174);
+  test_pow2 (174, 174, MPFR_RNDU, MPFR_RNDN, 174);
+  test_pow2 (176, 129, MPFR_RNDU, MPFR_RNDU, 174);
+  test_pow2 (176, 2, MPFR_RNDZ, MPFR_RNDZ, 174);
+  test_pow2 (176, 2, MPFR_RNDU, MPFR_RNDU, 176);
+
+  /* Tests for x = 2^i (E(x) = i+1) with error at most 1 = 2^0. */
+  for (n = 0; n < 100; n++)
+    {
+      i = (randlimb() % 200) + 4;
+      for (j = i - 2; j < i + 2; j++)
+        for (r1 = 0; r1 < MPFR_RND_MAX; r1++)
+          for (r2 = 0; r2 < MPFR_RND_MAX; r2++)
+            for (k = MPFR_PREC_MIN; k <= i + 2; k++)
+              test_pow2 (i, k, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j);
+    }
 
   mpfr_clear (x);