[src/fma.c] Fixed various bugs related to internal overflows/underflows.

[tests/tfma.c] Added tests.
(merged changesets r12393-12405,12583-12623 on these files from the
trunk)

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

2 files changed, 604 insertions, 226 deletions

diff --git a/src/fma.c b/src/fma.c
index 901fdb051..0706d4acb 100644
--- a/src/fma.c
+++ b/src/fma.c
@@ -225,194 +225,73 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
 
       if (MPFR_IS_INF (u))  /* overflow */
         {
+          int sign_u = MPFR_SIGN (u);
+
           MPFR_LOG_MSG (("Overflow on x*y\n", 0));
+          MPFR_GROUP_CLEAR (group);  /* we no longer need u */
 
           /* Let's eliminate the obvious case where x*y and z have the
              same sign. No possible cancellation -> real overflow.
              Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3,
-             then |x*y| >= 2^(emax+1), and |x*y + z| >= 2^emax. This case
+             then |x*y| >= 2^(emax+1), and |x*y + z| > 2^emax. This case
              is also an overflow. */
-          if (MPFR_SIGN (u) == MPFR_SIGN (z) || e >= __gmpfr_emax + 3)
+          if (sign_u == MPFR_SIGN (z) || e >= __gmpfr_emax + 3)
             {
-              MPFR_GROUP_CLEAR (group);
               MPFR_SAVE_EXPO_FREE (expo);
-              return mpfr_overflow (s, rnd_mode, MPFR_SIGN (z));
+              return mpfr_overflow (s, rnd_mode, sign_u);
             }
-
-          /* E(x) + E(y) <= emax+2, therefore |x*y| < 2^(emax+2), and
-             (x/4)*y does not overflow (let's recall that the result
-             is exact with an unbounded exponent range). It does not
-             underflow either, because x*y overflows and the exponent
-             range is large enough. */
-          inexact = mpfr_div_2ui (u, x, 2, MPFR_RNDN);
-          MPFR_ASSERTN (inexact == 0);
-          inexact = mpfr_mul (u, u, y, MPFR_RNDN);
-          MPFR_ASSERTN (inexact == 0);
-
-          /* Now, we need to add z/4... But it may underflow! */
-          {
-            mpfr_t zo4;
-            mpfr_srcptr zz;
-            MPFR_BLOCK_DECL (flags);
-
-            if (MPFR_GET_EXP (u) > MPFR_GET_EXP (z) &&
-                MPFR_GET_EXP (u) - MPFR_GET_EXP (z) > MPFR_PREC (u))
-              {
-                /* |z| < ulp(u)/2, therefore one can use z instead of z/4. */
-                zz = z;
-              }
-            else
-              {
-                mpfr_init2 (zo4, MPFR_PREC (z));
-                if (mpfr_div_2ui (zo4, z, 2, MPFR_RNDZ))
-                  {
-                    /* The division by 4 underflowed! */
-                    MPFR_ASSERTN (0); /* TODO... */
-                  }
-                zz = zo4;
-              }
-
-            /* Let's recall that u = x*y/4 and zz = z/4 (or z if the
-               following addition would give the same result). */
-            MPFR_BLOCK (flags, inexact = mpfr_add (s, u, zz, rnd_mode));
-            /* u and zz have different signs, so that an overflow
-               is not possible. But an underflow is theoretically
-               possible! */
-            if (MPFR_UNDERFLOW (flags))
-              {
-                MPFR_ASSERTN (zz != z);
-                MPFR_ASSERTN (0); /* TODO... */
-                mpfr_clears (zo4, u, (mpfr_ptr) 0);
-              }
-            else
-              {
-                int inex2;
-
-                if (zz != z)
-                  mpfr_clear (zo4);
-                MPFR_GROUP_CLEAR (group);
-                MPFR_ASSERTN (! MPFR_OVERFLOW (flags));
-                inex2 = mpfr_mul_2ui (s, s, 2, rnd_mode);
-                if (inex2)  /* overflow */
-                  {
-                    inexact = inex2;
-                    MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
-                  }
-                goto end;
-              }
-          }
         }
-      else  /* underflow: one has |xy| < 2^(emin-1). */
+      else  /* underflow: one has |x*y| < 2^(emin-1). */
         {
-          unsigned long scale = 0;
-          mpfr_t scaled_z;
-          mpfr_srcptr new_z;
-          mpfr_exp_t diffexp;
-          mpfr_prec_t pzs;
-          int xy_underflows;
-
           MPFR_LOG_MSG (("Underflow on x*y\n", 0));
 
-          /* Let's scale z so that ulp(z) > 2^emin and ulp(s) > 2^emin
-             (the + 1 on MPFR_PREC (s) is necessary because the exponent
-             of the result can be EXP(z) - 1). */
-          diffexp = MPFR_GET_EXP (z) - __gmpfr_emin;
-          pzs = MAX (MPFR_PREC (z), MPFR_PREC (s) + 1);
-          MPFR_LOG_MSG (("diffexp=%" MPFR_EXP_FSPEC "d pzs=%Pd\n",
-                         diffexp, pzs));
-          if (diffexp <= pzs)
-            {
-              mpfr_uexp_t uscale;
-              mpfr_t scaled_v;
-              MPFR_BLOCK_DECL (flags);
-
-              uscale = (mpfr_uexp_t) pzs - diffexp + 1;
-              MPFR_ASSERTN (uscale > 0);
-              MPFR_ASSERTN (uscale <= ULONG_MAX);
-              scale = uscale;
-              mpfr_init2 (scaled_z, MPFR_PREC (z));
-              inexact = mpfr_mul_2ui (scaled_z, z, scale, MPFR_RNDN);
-              MPFR_ASSERTN (inexact == 0);  /* TODO: overflow case */
-              new_z = scaled_z;
-              /* Now we need to recompute u = xy * 2^scale. */
-              MPFR_BLOCK (flags,
-                          if (MPFR_GET_EXP (x) < MPFR_GET_EXP (y))
-                            {
-                              mpfr_init2 (scaled_v, precx);
-                              mpfr_mul_2ui (scaled_v, x, scale, MPFR_RNDN);
-                              mpfr_mul (u, scaled_v, y, MPFR_RNDN);
-                            }
-                          else
-                            {
-                              mpfr_init2 (scaled_v, precy);
-                              mpfr_mul_2ui (scaled_v, y, scale, MPFR_RNDN);
-                              mpfr_mul (u, x, scaled_v, MPFR_RNDN);
-                            });
-              mpfr_clear (scaled_v);
-              MPFR_ASSERTN (! MPFR_OVERFLOW (flags));
-              xy_underflows = MPFR_UNDERFLOW (flags);
-            }
-          else
+          /* Easy cases: when 2^(emin-1) <= 1/2 * min(ulp(z),ulp(s)),
+             one can replace x*y by sign(x*y) * 2^(emin-1). Note that
+             this is even true in case of equality for MPFR_RNDN thanks
+             to the even-rounding rule.
+             The + 1 on MPFR_PREC (s) is necessary because the exponent
+             of the result can be EXP(z) - 1. */
+          if (MPFR_GET_EXP (z) - __gmpfr_emin >=
+              MAX (MPFR_PREC (z), MPFR_PREC (s) + 1))
             {
-              new_z = z;
-              xy_underflows = 1;
-            }
-
-          MPFR_LOG_MSG (("scale=%lu xy_underflows=%d\n",
-                         scale, xy_underflows));
-
-          if (xy_underflows)
-            {
-              /* Let's replace xy by sign(xy) * 2^(emin-1). */
               MPFR_PREC (u) = MPFR_PREC_MIN;
               mpfr_setmin (u, __gmpfr_emin);
               MPFR_SET_SIGN (u, MPFR_MULT_SIGN (MPFR_SIGN (x),
                                                 MPFR_SIGN (y)));
+              mpfr_clear_flags ();
+              goto add;
             }
 
-          {
-            MPFR_BLOCK_DECL (flags);
-
-            MPFR_BLOCK (flags, inexact = mpfr_add (s, u, new_z, rnd_mode));
-            MPFR_LOG_MSG (("inexact=%d\n", inexact));
-            MPFR_GROUP_CLEAR (group);
-            if (scale != 0)
-              {
-                int inex2;
-
-                mpfr_clear (scaled_z);
-                /* Here an overflow is theoretically possible, in which case
-                   the result may be wrong, hence the assert. An underflow
-                   is not possible, but let's check that anyway. */
-                MPFR_ASSERTN (! MPFR_OVERFLOW (flags));  /* TODO... */
-                MPFR_ASSERTN (! MPFR_UNDERFLOW (flags));  /* not possible */
-                if (rnd_mode == MPFR_RNDN &&
-                    MPFR_GET_EXP (s) == __gmpfr_emin - 1 + scale &&
-                    mpfr_powerof2_raw (s))
-                  {
-                    MPFR_LOG_MSG (("Double rounding\n", 0));
-                    rnd_mode = (MPFR_IS_NEG (s) ? inexact <= 0 : inexact >= 0)
-                      ? MPFR_RNDZ : MPFR_RNDA;
-                  }
-                inex2 = mpfr_div_2ui (s, s, scale, rnd_mode);
-                MPFR_LOG_MSG (("inex2=%d\n", inex2));
-                if (inex2)  /* underflow */
-                  inexact = inex2;
-              }
-          }
-
-          /* FIXME/TODO: I'm not sure that the following is correct.
-             Check for possible spurious exceptions due to intermediate
-             computations. */
-          MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
-          goto end;
+          MPFR_GROUP_CLEAR (group);  /* we no longer need u */
         }
+
+      /* Let's use UBF to resolve the overflow/underflow issues. */
+      {
+        mpfr_ubf_t uu;
+        mp_size_t un;
+        mpfr_limb_ptr up;
+        MPFR_TMP_DECL(marker);
+
+        MPFR_LOG_MSG (("Use UBF\n", 0));
+
+        MPFR_TMP_MARK (marker);
+        un = MPFR_LIMB_SIZE (x) + MPFR_LIMB_SIZE (y);
+        MPFR_TMP_INIT (up, uu, (mpfr_prec_t) un * GMP_NUMB_BITS, un);
+        mpfr_ubf_mul_exact (uu, x, y);
+        mpfr_clear_flags ();
+        inexact = mpfr_add (s, (mpfr_srcptr) uu, z, rnd_mode);
+        MPFR_UBF_CLEAR_EXP (uu);
+        MPFR_TMP_FREE (marker);
+      }
+    }
+  else
+    {
+    add:
+      inexact = mpfr_add (s, u, z, rnd_mode);
+      MPFR_GROUP_CLEAR (group);
     }
 
-  inexact = mpfr_add (s, u, z, rnd_mode);
-  MPFR_GROUP_CLEAR (group);
   MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
- end:
   MPFR_SAVE_EXPO_FREE (expo);
   return mpfr_check_range (s, inexact, rnd_mode);
 }
diff --git a/tests/tfma.c b/tests/tfma.c
index f3a3b30e7..c09986b0a 100644
--- a/tests/tfma.c
+++ b/tests/tfma.c
@@ -196,6 +196,238 @@ test_overflow2 (void)
 }
 
 static void
+test_overflow3 (void)
+{
+  mpfr_t x, y, z, r;
+  int inex;
+  mpfr_prec_t p = 8;
+  mpfr_flags_t ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT, flags;
+  int i, j, k;
+  unsigned int neg;
+
+  mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0);
+  for (i = 0; i < 2; i++)
+    {
+      mpfr_init2 (r, 2 * p + i);
+      mpfr_set_ui_2exp (x, 1, mpfr_get_emax () - 1, MPFR_RNDN);
+      mpfr_set_ui (y, 2, MPFR_RNDN);       /* y = 2 */
+      for (j = 1; j <= 2; j++)
+        for (k = 0; k <= 1; k++)
+          {
+            mpfr_set_si_2exp (z, -1, mpfr_get_emax () - mpfr_get_prec (r) - j,
+                              MPFR_RNDN);
+            if (k)
+              mpfr_nextabove (z);
+            for (neg = 0; neg <= 3; neg++)
+              {
+                mpfr_clear_flags ();
+                /* (The following applies for neg = 0 or 2, all the signs
+                   need to be reversed for neg = 1 or 3.)
+                   We have x*y = 2^emax and
+                   z = - 2^(emax-2p-i-j) * (1-k*2^(-p)), thus
+                   x*y+z = 2^emax - 2^(emax-2p-i-j) + k*2^(emax-3p-i-j)
+                   should overflow. Indeed it is >= the midpoint of
+                   2^emax - 2^(emax-2p-i) and 2^emax, the midpoint
+                   being obtained for j = 1 and k = 0. */
+                inex = mpfr_fma (r, x, y, z, MPFR_RNDN);
+                flags = __gmpfr_flags;
+                if (! mpfr_inf_p (r) || flags != ex_flags ||
+                    ((neg & 1) == 0 ?
+                     (inex <= 0 || MPFR_IS_NEG (r)) :
+                     (inex >= 0 || MPFR_IS_POS (r))))
+                  {
+                    printf ("Error in test_overflow3 for "
+                            "i=%d j=%d k=%d neg=%u\n", i, j, k, neg);
+                    printf ("Expected %c@Inf@\n  with inex %c 0 and flags:",
+                            (neg & 1) == 0 ? '+' : '-',
+                            (neg & 1) == 0 ? '>' : '<');
+                    flags_out (ex_flags);
+                    printf ("Got      ");
+                    mpfr_dump (r);
+                    printf ("  with inex = %d and flags:", inex);
+                    flags_out (flags);
+                    exit (1);
+                  }
+                if (neg == 0 || neg == 2)
+                  mpfr_neg (x, x, MPFR_RNDN);
+                if (neg == 1 || neg == 3)
+                  mpfr_neg (y, y, MPFR_RNDN);
+                mpfr_neg (z, z, MPFR_RNDN);
+              }  /* neg */
+          }  /* k */
+      mpfr_clear (r);
+    }  /* i */
+  mpfr_clears (x, y, z, (mpfr_ptr) 0);
+}
+
+static void
+test_overflow4 (void)
+{
+  mpfr_t x, y, z, r1, r2;
+  mpfr_exp_t emax, e;
+  mpfr_prec_t px;
+  mpfr_flags_t flags1, flags2;
+  int inex1, inex2;
+  int ei, i, j;
+  int below;
+  unsigned int neg;
+
+  emax = mpfr_get_emax ();
+
+  mpfr_init2 (y, MPFR_PREC_MIN);
+  mpfr_set_ui (y, 2, MPFR_RNDN);  /* y = 2 */
+
+  mpfr_init2 (z, 8);
+
+  for (px = 17; px < 256; px *= 2)
+    {
+      mpfr_init2 (x, px);
+      mpfr_inits2 (px - 8, r1, r2, (mpfr_ptr) 0);
+      for (ei = 0; ei <= 1; ei++)
+        {
+          e = ei ? emax : 0;
+          mpfr_set_ui_2exp (x, 1, e - 1, MPFR_RNDN);
+          mpfr_nextabove (x);  /* x = 2^(e - 1) + 2^(e - px) */
+          /* x*y = 2^e + 2^(e - px + 1), which internally overflows
+             when e = emax. */
+          for (i = -4; i <= 4; i++)
+            for (j = 2; j <= 3; j++)
+              {
+                mpfr_set_si_2exp (z, -j, e - px + i, MPFR_RNDN);
+                /* If |z| <= 2^(e - px + 1), then x*y + z >= 2^e and
+                   RZ(x*y + z) = 2^e with an unbounded exponent range.
+                   If |z| > 2^(e - px + 1), then RZ(x*y + z) is the
+                   predecessor of 2^e (since |z| < ulp(r)/2); this
+                   occurs when i > 0 and when i = 0 and j > 2 */
+                mpfr_set_ui_2exp (r1, 1, e - 1, MPFR_RNDN);
+                below = i > 0 || (i == 0 && j > 2);
+                if (below)
+                  mpfr_nextbelow (r1);
+                mpfr_clear_flags ();
+                inex1 = mpfr_mul_2ui (r1, r1, 1, MPFR_RNDZ);
+                if (below || e < emax)
+                  {
+                    inex1 = i == 0 && j == 2 ? 0 : -1;
+                    flags1 = inex1 ? MPFR_FLAGS_INEXACT : 0;
+                  }
+                else
+                  {
+                    MPFR_ASSERTN (inex1 < 0);
+                    flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
+                    MPFR_ASSERTN (flags1 == __gmpfr_flags);
+                  }
+                for (neg = 0; neg <= 3; neg++)
+                  {
+                    mpfr_clear_flags ();
+                    inex2 = mpfr_fma (r2, x, y, z, MPFR_RNDZ);
+                    flags2 = __gmpfr_flags;
+                    if (! (mpfr_equal_p (r1, r2) &&
+                           SAME_SIGN (inex1, inex2) &&
+                           flags1 == flags2))
+                      {
+                        printf ("Error in test_overflow4 for "
+                                "px=%d ei=%d i=%d j=%d neg=%u\n",
+                                (int) px, ei, i, j, neg);
+                        printf ("Expected ");
+                        mpfr_dump (r1);
+                        printf ("with inex = %d and flags:", inex1);
+                        flags_out (flags1);
+                        printf ("Got      ");
+                        mpfr_dump (r2);
+                        printf ("with inex = %d and flags:", inex2);
+                        flags_out (flags2);
+                        exit (1);
+                      }
+                    if (neg == 0 || neg == 2)
+                      mpfr_neg (x, x, MPFR_RNDN);
+                    if (neg == 1 || neg == 3)
+                      mpfr_neg (y, y, MPFR_RNDN);
+                    mpfr_neg (z, z, MPFR_RNDN);
+                    mpfr_neg (r1, r1, MPFR_RNDN);
+                    inex1 = - inex1;
+                  }
+              }
+        }
+      mpfr_clears (x, r1, r2, (mpfr_ptr) 0);
+    }
+
+  mpfr_clears (y, z, (mpfr_ptr) 0);
+}
+
+static void
+test_overflow5 (void)
+{
+  mpfr_t x, y, z, r1, r2;
+  mpfr_exp_t emax;
+  int inex1, inex2;
+  int i, rnd;
+  unsigned int neg, negr;
+
+  emax = mpfr_get_emax ();
+
+  mpfr_init2 (x, 123);
+  mpfr_init2 (y, 45);
+  mpfr_init2 (z, 67);
+  mpfr_inits2 (89, r1, r2, (mpfr_ptr) 0);
+
+  mpfr_set_ui_2exp (x, 1, emax - 1, MPFR_RNDN);
+
+  for (i = 3; i <= 17; i++)
+    {
+      mpfr_set_ui (y, i, MPFR_RNDN);
+      mpfr_set_ui_2exp (z, 1, emax - 1, MPFR_RNDN);
+      for (neg = 0; neg < 8; neg++)
+        {
+          mpfr_setsign (x, x, neg & 1, MPFR_RNDN);
+          mpfr_setsign (y, y, neg & 2, MPFR_RNDN);
+          mpfr_setsign (z, z, neg & 4, MPFR_RNDN);
+
+          /* |x*y + z| = (i +/- 1) * 2^(emax - 1) >= 2^emax (overflow)
+             and x*y + z has the same sign as x*y. */
+          negr = (neg ^ (neg >> 1)) & 1;
+
+          RND_LOOP (rnd)
+            {
+              mpfr_set_inf (r1, 1);
+              if (MPFR_IS_LIKE_RNDZ ((mpfr_rnd_t) rnd, negr))
+                {
+                  mpfr_nextbelow (r1);
+                  inex1 = -1;
+                }
+              else
+                inex1 = 1;
+
+              if (negr)
+                {
+                  mpfr_neg (r1, r1, MPFR_RNDN);
+                  inex1 = - inex1;
+                }
+
+              mpfr_clear_flags ();
+              inex2 = mpfr_fma (r2, x, y, z, (mpfr_rnd_t) rnd);
+              MPFR_ASSERTN (__gmpfr_flags ==
+                            (MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW));
+
+              if (! (mpfr_equal_p (r1, r2) && SAME_SIGN (inex1, inex2)))
+                {
+                  printf ("Error in test_overflow5 for i=%d neg=%u %s\n",
+                          i, neg, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+                  printf ("Expected ");
+                  mpfr_dump (r1);
+                  printf ("with inex = %d\n", inex1);
+                  printf ("Got      ");
+                  mpfr_dump (r2);
+                  printf ("with inex = %d\n", inex2);
+                  exit (1);
+                }
+            }  /* rnd */
+        }  /* neg */
+    }  /* i */
+
+  mpfr_clears (x, y, z, r1, r2, (mpfr_ptr) 0);
+}
+
+static void
 test_underflow1 (void)
 {
   mpfr_t x, y, z, r;
@@ -281,59 +513,128 @@ test_underflow1 (void)
 static void
 test_underflow2 (void)
 {
-  mpfr_t x, y, z, r;
-  int b, i, inex, same, err = 0;
+  mpfr_t x, y, z, r1, r2;
+  int e, b, i, prec = 32, pz, inex;
+  unsigned int neg;
 
-  mpfr_inits2 (32, x, y, z, r, (mpfr_ptr) 0);
+  mpfr_init2 (x, MPFR_PREC_MIN);
+  mpfr_inits2 (prec, y, z, r1, r2, (mpfr_ptr) 0);
 
-  mpfr_set_si_2exp (z, 1, mpfr_get_emin (), MPFR_RNDN);   /* z = 2^emin */
-  mpfr_set_si_2exp (x, 1, mpfr_get_emin (), MPFR_RNDN);   /* x = 2^emin */
+  mpfr_set_si_2exp (x, 1, mpfr_get_emin () - 1, MPFR_RNDN);
+  /* x = 2^(emin-1) */
 
-  for (b = 0; b <= 1; b++)
+  for (e = -1; e <= prec + 2; e++)
     {
-      for (i = 15; i <= 17; i++)
+      mpfr_set (z, x, MPFR_RNDN);
+      /* z = x = 2^(emin+e) */
+      for (b = 0; b <= 1; b++)
         {
-          mpfr_set_si_2exp (y, i, -4 - MPFR_PREC (z), MPFR_RNDN);
-          /*  z = 1.000...00b
-           * xy =            01111
-           *   or            10000
-           *   or            10001
-           */
-          mpfr_clear_flags ();
-          inex = mpfr_fma (r, x, y, z, MPFR_RNDN);
-#define STRTU2 "Error in test_underflow2 (b = %d, i = %d)\n  "
-          if (__gmpfr_flags != MPFR_FLAGS_INEXACT)
+          for (pz = prec - 4 * (b == 0); pz <= prec + 4; pz++)
             {
-              printf (STRTU2 "flags = %u instead of %u\n", b, i,
-                      (unsigned int) __gmpfr_flags,
-                      (unsigned int) MPFR_FLAGS_INEXACT);
-              err = 1;
-            }
-          same = i == 15 || (i == 16 && b == 0);
-          if (same ? (inex >= 0) : (inex <= 0))
-            {
-              printf (STRTU2 "incorrect ternary value (%d instead of %c 0)\n",
-                      b, i, inex, same ? '<' : '>');
-              err = 1;
-            }
-          mpfr_set (y, z, MPFR_RNDN);
-          if (!same)
-            mpfr_nextabove (y);
-          if (! mpfr_equal_p (r, y))
-            {
-              printf (STRTU2 "expected ", b, i);
-              mpfr_dump (y);
-              printf ("  got      ");
-              mpfr_dump (r);
-              err = 1;
-            }
-        }
-      mpfr_nextabove (z);
-    }
-
-  if (err)
-    exit (1);
-  mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
+              inex = mpfr_prec_round (z, pz, MPFR_RNDZ);
+              MPFR_ASSERTN (inex == 0);
+              for (i = 15; i <= 17; i++)
+                {
+                  mpfr_flags_t flags1, flags2;
+                  int inex1, inex2;
+
+                  mpfr_set_si_2exp (y, i, -4 - prec, MPFR_RNDN);
+
+                  /*      <--- r --->
+                   *  z = 1.000...00b   with b = 0 or 1
+                   * xy =            01111  (i = 15)
+                   *   or            10000  (i = 16)
+                   *   or            10001  (i = 17)
+                   *
+                   * x, y, and z will be modified to test the different sign
+                   * combinations. In the case b = 0 (i.e. |z| is a power of
+                   * 2) and x*y has a different sign from z, then y will be
+                   * divided by 2, so that i = 16 corresponds to a midpoint.
+                   */
+
+                  for (neg = 0; neg < 8; neg++)
+                    {
+                      int xyneg, prev_binade;
+
+                      mpfr_setsign (x, x, neg & 1, MPFR_RNDN);
+                      mpfr_setsign (y, y, neg & 2, MPFR_RNDN);
+                      mpfr_setsign (z, z, neg & 4, MPFR_RNDN);
+
+                      xyneg = (neg ^ (neg >> 1)) & 1;  /* true iff x*y < 0 */
+
+                      /* If a change of binade occurs by adding x*y to z
+                         exactly, then take into account the fact that the
+                         midpoint has a different exponent. */
+                      prev_binade = b == 0 && (xyneg ^ MPFR_IS_NEG (z));
+                      if (prev_binade)
+                        mpfr_div_2ui (y, y, 1, MPFR_RNDN);
+
+                      mpfr_set (r1, z, MPFR_RNDN);
+                      flags1 = MPFR_FLAGS_INEXACT;
+
+                      if (e == -1 && i == 17 && b == 0 &&
+                          (xyneg ^ (neg >> 2)) != 0)
+                        {
+                          /* Special underflow case. */
+                          flags1 |= MPFR_FLAGS_UNDERFLOW;
+                          inex1 = xyneg ? 1 : -1;
+                        }
+                      else if (i == 15 || (i == 16 && b == 0))
+                        {
+                          /* round toward z */
+                          inex1 = xyneg ? 1 : -1;
+                        }
+                      else if (xyneg)
+                        {
+                          /* round away from z, with x*y < 0 */
+                          mpfr_nextbelow (r1);
+                          inex1 = -1;
+                        }
+                      else
+                        {
+                          /* round away from z, with x*y > 0 */
+                          mpfr_nextabove (r1);
+                          inex1 = 1;
+                        }
+
+                      mpfr_clear_flags ();
+                      inex2 = mpfr_fma (r2, x, y, z, MPFR_RNDN);
+                      flags2 = __gmpfr_flags;
+
+                      if (! (mpfr_equal_p (r1, r2) &&
+                             SAME_SIGN (inex1, inex2) &&
+                             flags1 == flags2))
+                        {
+                          printf ("Error in test_underflow2 for "
+                                  "e=%d b=%d pz=%d i=%d neg=%u\n",
+                                  e, b, pz, i, neg);
+                          printf ("Expected ");
+                          mpfr_dump (r1);
+                          printf ("with inex = %d and flags:", inex1);
+                          flags_out (flags1);
+                          printf ("Got      ");
+                          mpfr_dump (r2);
+                          printf ("with inex = %d and flags:", inex2);
+                          flags_out (flags2);
+                          exit (1);
+                        }
+
+                      /* Restore y. */
+                      if (prev_binade)
+                        mpfr_mul_2ui (y, y, 1, MPFR_RNDN);
+                    }  /* neg */
+                }  /* i */
+            }  /* pz */
+
+          inex = mpfr_prec_round (z, prec, MPFR_RNDZ);
+          MPFR_ASSERTN (inex == 0);
+          MPFR_SET_POS (z);
+          mpfr_nextabove (z);
+        }  /* b */
+      mpfr_mul_2ui (x, x, 1, MPFR_RNDN);
+    }  /* e */
+
+  mpfr_clears (x, y, z, r1, r2, (mpfr_ptr) 0);
 }
 
 static void
@@ -397,6 +698,185 @@ test_underflow3 (int n)
   mpfr_clears (x, y, z, t1, t2, (mpfr_ptr) 0);
 }
 
+/* Test s = x*y + z with PREC(z) > PREC(s) + 1, x*y underflows, where
+   z + x*y and z + sign(x*y) * 2^(emin-1) do not give the same result.
+     x = 2^emin
+     y = 2^(-8)
+     z = 2^emin * (2^PREC(s) + k - 2^(-1))
+   with k = 3 for MPFR_RNDN and k = 2 for the directed rounding modes.
+   Also test the opposite versions with neg != 0.
+*/
+static void
+test_underflow4 (void)
+{
+  mpfr_t x, y, z, s1, s2;
+  mpfr_prec_t ps = 32;
+  int inex, rnd;
+
+  mpfr_inits2 (MPFR_PREC_MIN, x, y, (mpfr_ptr) 0);
+  mpfr_inits2 (ps, s1, s2, (mpfr_ptr) 0);
+  mpfr_init2 (z, ps + 2);
+
+  inex = mpfr_set_si_2exp (x, 1, mpfr_get_emin (), MPFR_RNDN);
+  MPFR_ASSERTN (inex == 0);
+  inex = mpfr_set_si_2exp (y, 1, -8, MPFR_RNDN);
+  MPFR_ASSERTN (inex == 0);
+
+  RND_LOOP_NO_RNDF (rnd)
+    {
+      mpfr_flags_t flags1, flags2;
+      int inex1, inex2;
+      unsigned int neg;
+
+      inex = mpfr_set_si_2exp (z, 1 << 1, ps, MPFR_RNDN);
+      MPFR_ASSERTN (inex == 0);
+      inex = mpfr_sub_ui (z, z, 1, MPFR_RNDN);
+      MPFR_ASSERTN (inex == 0);
+      inex = mpfr_div_2ui (z, z, 1, MPFR_RNDN);
+      MPFR_ASSERTN (inex == 0);
+      inex = mpfr_add_ui (z, z, rnd == MPFR_RNDN ? 3 : 2, MPFR_RNDN);
+      MPFR_ASSERTN (inex == 0);
+      inex = mpfr_mul (z, z, x, MPFR_RNDN);
+      MPFR_ASSERTN (inex == 0);
+
+      for (neg = 0; neg <= 3; neg++)
+        {
+          inex1 = mpfr_set (s1, z, (mpfr_rnd_t) rnd);
+          flags1 = MPFR_FLAGS_INEXACT;
+
+          mpfr_clear_flags ();
+          inex2 = mpfr_fma (s2, x, y, z, (mpfr_rnd_t) rnd);
+          flags2 = __gmpfr_flags;
+
+          if (! (mpfr_equal_p (s1, s2) &&
+                 SAME_SIGN (inex1, inex2) &&
+                 flags1 == flags2))
+            {
+              printf ("Error in test_underflow4 for neg=%u %s\n",
+                      neg, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+              printf ("Expected ");
+              mpfr_dump (s1);
+              printf ("  with inex ~ %d, flags =", inex1);
+              flags_out (flags1);
+              printf ("Got      ");
+              mpfr_dump (s2);
+              printf ("  with inex ~ %d, flags =", inex2);
+              flags_out (flags2);
+              exit (1);
+            }
+
+          if (neg == 0 || neg == 2)
+            mpfr_neg (x, x, MPFR_RNDN);
+          if (neg == 1 || neg == 3)
+            mpfr_neg (y, y, MPFR_RNDN);
+          mpfr_neg (z, z, MPFR_RNDN);
+        }
+    }
+
+  mpfr_clears (x, y, z, s1, s2, (mpfr_ptr) 0);
+}
+
+/* Test s = x*y + z on:
+     x = +/- 2^emin
+     y = +/- 2^(-3)
+     z = +/- 2^(emin + PREC(s)) and MPFR numbers close to this value.
+   with PREC(z) from PREC(s) - 2 to PREC(s) + 8.
+*/
+static void
+test_underflow5 (void)
+{
+  mpfr_t w, x, y, z, s1, s2, t;
+  mpfr_exp_t emin;
+  mpfr_prec_t ps = 32;
+  int inex, i, j, rnd;
+  unsigned int neg;
+
+  mpfr_inits2 (MPFR_PREC_MIN, w, x, y, (mpfr_ptr) 0);
+  mpfr_inits2 (ps, s1, s2, (mpfr_ptr) 0);
+  mpfr_init2 (t, ps + 9);
+
+  emin = mpfr_get_emin ();
+
+  inex = mpfr_set_si_2exp (w, 1, emin, MPFR_RNDN);  /* w = 2^emin */
+  MPFR_ASSERTN (inex == 0);
+  inex = mpfr_set_si (x, 1, MPFR_RNDN);
+  MPFR_ASSERTN (inex == 0);
+  inex = mpfr_set_si_2exp (y, 1, -3, MPFR_RNDN);
+  MPFR_ASSERTN (inex == 0);
+
+  for (i = -2; i <= 8; i++)
+    {
+      mpfr_init2 (z, ps + i);
+      inex = mpfr_set_si_2exp (z, 1, ps, MPFR_RNDN);
+      MPFR_ASSERTN (inex == 0);
+
+      for (j = 1; j <= 32; j++)
+        mpfr_nextbelow (z);
+
+      for (j = -32; j <= 32; j++)
+        {
+          for (neg = 0; neg < 8; neg++)
+            {
+              mpfr_setsign (x, x, neg & 1, MPFR_RNDN);
+              mpfr_setsign (y, y, neg & 2, MPFR_RNDN);
+              mpfr_setsign (z, z, neg & 4, MPFR_RNDN);
+
+              inex = mpfr_fma (t, x, y, z, MPFR_RNDN);
+              MPFR_ASSERTN (inex == 0);
+
+              inex = mpfr_mul (x, x, w, MPFR_RNDN);
+              MPFR_ASSERTN (inex == 0);
+              inex = mpfr_mul (z, z, w, MPFR_RNDN);
+              MPFR_ASSERTN (inex == 0);
+
+              RND_LOOP_NO_RNDF (rnd)
+                {
+                  mpfr_flags_t flags1, flags2;
+                  int inex1, inex2;
+
+                  mpfr_clear_flags ();
+                  inex1 = mpfr_mul (s1, w, t, (mpfr_rnd_t) rnd);
+                  flags1 = __gmpfr_flags;
+
+                  mpfr_clear_flags ();
+                  inex2 = mpfr_fma (s2, x, y, z, (mpfr_rnd_t) rnd);
+                  flags2 = __gmpfr_flags;
+
+                  if (! (mpfr_equal_p (s1, s2) &&
+                         SAME_SIGN (inex1, inex2) &&
+                         flags1 == flags2))
+                    {
+                      printf ("Error in test_underflow5 on "
+                              "i=%d j=%d neg=%u %s\n", i, j, neg,
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+                      printf ("Expected ");
+                      mpfr_dump (s1);
+                      printf ("  with inex ~ %d, flags =", inex1);
+                      flags_out (flags1);
+                      printf ("Got      ");
+                      mpfr_dump (s2);
+                      printf ("  with inex ~ %d, flags =", inex2);
+                      flags_out (flags2);
+                      exit (1);
+                    }
+                }  /* rnd */
+
+              inex = mpfr_div (x, x, w, MPFR_RNDN);
+              MPFR_ASSERTN (inex == 0);
+              inex = mpfr_div (z, z, w, MPFR_RNDN);
+              MPFR_ASSERTN (inex == 0);
+            }  /* neg */
+
+          MPFR_SET_POS (z);  /* restore the value before the loop on neg */
+          mpfr_nextabove (z);
+        }  /* j */
+
+      mpfr_clear (z);
+    }  /* i */
+
+  mpfr_clears (w, x, y, s1, s2, t, (mpfr_ptr) 0);
+}
+
 static void
 bug20101018 (void)
 {
@@ -533,6 +1013,7 @@ main (int argc, char *argv[])
 {
   mpfr_t x, y, z, s;
   mpfr_exp_t emin, emax;
+  int i;
 
   tests_start_mpfr ();
 
@@ -823,21 +1304,39 @@ main (int argc, char *argv[])
 
   test_exact ();
 
-  test_overflow1 ();
-  test_overflow2 ();
-  test_underflow1 ();
-  test_underflow2 ();
-  test_underflow3 (1);
-
-  set_emin (MPFR_EMIN_MIN);
-  set_emax (MPFR_EMAX_MAX);
-  test_overflow1 ();
-  test_overflow2 ();
-  test_underflow1 ();
-  test_underflow2 ();
-  test_underflow3 (2);
-  set_emin (emin);
-  set_emax (emax);
+  for (i = 0; i <= 2; i++)
+    {
+      if (i == 0)
+        {
+          /* corresponds to the generic code + mpfr_check_range */
+          set_emin (-1024);
+          set_emax (1024);
+        }
+      else if (i == 1)
+        {
+          set_emin (MPFR_EMIN_MIN);
+          set_emax (MPFR_EMAX_MAX);
+        }
+      else
+        {
+          MPFR_ASSERTN (i == 2);
+          if (emin == MPFR_EMIN_MIN && emax == MPFR_EMAX_MAX)
+            break;
+          set_emin (emin);
+          set_emax (emax);
+        }
+
+      test_overflow1 ();
+      test_overflow2 ();
+      test_overflow3 ();
+      test_overflow4 ();
+      test_overflow5 ();
+      test_underflow1 ();
+      test_underflow2 ();
+      test_underflow3 (i);
+      test_underflow4 ();
+      test_underflow5 ();
+    }
 
   tests_end_mpfr ();
   return 0;