Fixed double-rounding problem in mpfr_mul_2si, mpfr_div_2si

and mpfr_div_2ui rounded to nearest when the rounded result
in unbounded exponent range is 2^(emin - 2): the returned
result was zero instead of non-zero. Added testcases.
[Merged changesets 5503 5510 5511 5512 5517 from the trunk]

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

4 files changed, 93 insertions, 3 deletions

diff --git a/div_2si.c b/div_2si.c
index 4ead0a9c7..b01e69c34 100644
--- a/div_2si.c
+++ b/div_2si.c
@@ -37,7 +37,8 @@ mpfr_div_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
         {
           if (rnd_mode == GMP_RNDN &&
               (__gmpfr_emin > MPFR_EMAX_MAX - (n - 1) ||
-               exp < __gmpfr_emin + (n - 1) || mpfr_powerof2_raw (y)))
+               exp < __gmpfr_emin + (n - 1) ||
+               (inexact >= 0 && mpfr_powerof2_raw (y))))
             rnd_mode = GMP_RNDZ;
           return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
         }
diff --git a/div_2ui.c b/div_2ui.c
index af19fdaec..bca99cd9f 100644
--- a/div_2ui.c
+++ b/div_2ui.c
@@ -55,7 +55,7 @@ mpfr_div_2ui (mpfr_ptr y, mpfr_srcptr x, unsigned long n, mp_rnd_t rnd_mode)
             if (rnd_mode == GMP_RNDN &&
                 (__gmpfr_emin > MPFR_EMAX_MAX - (long) (n - 1) ||
                  exp < __gmpfr_emin + (long) (n - 1) ||
-                 mpfr_powerof2_raw (y)))
+                 (inexact >= 0 && mpfr_powerof2_raw (y))))
               rnd_mode = GMP_RNDZ;
             return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
           }
diff --git a/mul_2si.c b/mul_2si.c
index 7a564b02b..d4d83893a 100644
--- a/mul_2si.c
+++ b/mul_2si.c
@@ -41,7 +41,8 @@ mpfr_mul_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
         {
           if (rnd_mode == GMP_RNDN &&
               (__gmpfr_emin > MPFR_EMAX_MAX + (n + 1) ||
-               exp < __gmpfr_emin - (n + 1) || mpfr_powerof2_raw (y)))
+               exp < __gmpfr_emin - (n + 1) ||
+               (inexact >= 0 && mpfr_powerof2_raw (y))))
             rnd_mode = GMP_RNDZ;
           return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
         }
diff --git a/tests/tmul_2exp.c b/tests/tmul_2exp.c
index 4dcca9455..3b1c23596 100644
--- a/tests/tmul_2exp.c
+++ b/tests/tmul_2exp.c
@@ -45,6 +45,92 @@ test_mul (int i, int div, mpfr_ptr y, mpfr_srcptr x,
     (exit (1), 0);
 }
 
+static void
+underflow (mp_exp_t e)
+{
+  mpfr_t x, y, z1, z2;
+  mp_exp_t emin;
+  int i, k;
+  int prec;
+  int rnd;
+  int div;
+  int inex1, inex2;
+  unsigned int flags1, flags2;
+
+  /* Test mul_2si(x, e - k), div_2si(x, k - e) and div_2ui(x, k - e)
+   * with emin = e, x = 1 + i/16, i in { -1, 0, 1 }, and k = 1 to 4,
+   * by comparing the result with the one of a simple division.
+   */
+  emin = mpfr_get_emin ();
+  set_emin (e);
+  mpfr_inits2 (8, x, y, (mpfr_ptr) 0);
+  for (i = 15; i <= 17; i++)
+    {
+      inex1 = mpfr_set_ui_2exp (x, i, -4, GMP_RNDN);
+      MPFR_ASSERTN (inex1 == 0);
+      for (prec = 6; prec >= 3; prec -= 3)
+        {
+          mpfr_inits2 (prec, z1, z2, (mpfr_ptr) 0);
+          RND_LOOP (rnd)
+            for (k = 1; k <= 4; k++)
+              {
+                /* The following one is assumed to be correct. */
+                inex1 = mpfr_mul_2si (y, x, e, GMP_RNDN);
+                MPFR_ASSERTN (inex1 == 0);
+                inex1 = mpfr_set_ui (z1, 1 << k, GMP_RNDN);
+                MPFR_ASSERTN (inex1 == 0);
+                mpfr_clear_flags ();
+                /* Do not use mpfr_div_ui to avoid the optimization
+                   by mpfr_div_2si. */
+                inex1 = mpfr_div (z1, y, z1, rnd);
+                flags1 = __gmpfr_flags;
+
+              for (div = 0; div <= 2; div++)
+                {
+                  mpfr_clear_flags ();
+                  inex2 = div == 0 ?
+                    mpfr_mul_2si (z2, x, e - k, rnd) : div == 1 ?
+                    mpfr_div_2si (z2, x, k - e, rnd) :
+                    mpfr_div_2ui (z2, x, k - e, rnd);
+                  flags2 = __gmpfr_flags;
+                  if (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                      mpfr_equal_p (z1, z2))
+                    continue;
+                  printf ("Error in underflow(");
+                  if (e == MPFR_EMIN_MIN)
+                    printf ("MPFR_EMIN_MIN");
+                  else if (e == emin)
+                    printf ("default emin");
+                  else
+                    printf ("%ld", e);
+                  printf (")\nwith %s, x = %d/16, prec = %d, k = %d, mpfr_%s",
+                          div == 0 ? "mul_2si" : div == 1 ? "div_2si" :
+                          "div_2ui", i, prec, k, mpfr_print_rnd_mode (rnd));
+                  printf ("\nExpected ");
+                  mpfr_out_str (stdout, 16, 0, z1, GMP_RNDN);
+                  printf (", inex = %d, flags = %u\n", SIGN (inex1), flags1);
+                  printf ("Got      ");
+                  mpfr_out_str (stdout, 16, 0, z2, GMP_RNDN);
+                  printf (", inex = %d, flags = %u\n", SIGN (inex2), flags2);
+                  exit (1);
+                }  /* div */
+              }  /* k */
+          mpfr_clears (z1, z2, (mpfr_ptr) 0);
+        }  /* prec */
+    }  /* i */
+  mpfr_clears (x, y, (mpfr_ptr) 0);
+  set_emin (emin);
+}
+
+static void
+underflow0 (void)
+{
+  underflow (-256);
+  if (mpfr_get_emin () != MPFR_EMIN_MIN)
+    underflow (mpfr_get_emin ());
+  underflow (MPFR_EMIN_MIN);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -134,6 +220,8 @@ main (int argc, char *argv[])
 
   mpfr_clears (w, z, (mpfr_ptr) 0);
 
+  underflow0 ();
+
   tests_end_mpfr ();
   return 0;
 }