[src/sub1.c] Fixed bug in mpfr_sub1 (real subtraction b - c, |b| > |c|):

  In MPFR_RNDN (rounding to nearest), when |b| is the midpoint between
  the maximum number and 2^emax (the maximum number + 1 ulp) and c is
  small, the obtained result is an infinity (with overflow) instead of
  ± the maximum number (no overflow). The cause is that an overflow is
  generated too early (in the rounding code).
[tests/tsub.c] Added test cases.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10383 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 139 insertions, 5 deletions

diff --git a/src/sub1.c b/src/sub1.c
index a0bee3d6a..1929bb4fc 100644
--- a/src/sub1.c
+++ b/src/sub1.c
@@ -106,16 +106,15 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
       /* A = S*ABS(B) +/- ulp(a) */
       MPFR_SET_EXP (a, MPFR_GET_EXP (b));
       MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq,
-                        rnd_mode, MPFR_SIGN (a),
-                        if (MPFR_UNLIKELY (++ MPFR_EXP (a) > __gmpfr_emax))
-                        inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)));
-      /* inexact = mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (a));  */
+                        rnd_mode, MPFR_SIGN (a), ++ MPFR_EXP (a));
       if (inexact == 0)
         {
           /* a = b (Exact)
              But we know it isn't (Since we have to remove `c')
              So if we round to Zero, we have to remove one ulp.
              Otherwise the result is correctly rounded. */
+          /* An overflow is not possible. */
+          MPFR_ASSERTD (MPFR_EXP (a) <= __gmpfr_emax);
           if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
             {
               mpfr_nexttozero (a);
@@ -146,9 +145,14 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
              i.e. inexact= MPFR_EVEN_INEX */
           if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX * MPFR_INT_SIGN (a)))
             {
-              mpfr_nexttozero (a);
+              if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+                mpfr_setmax (a, __gmpfr_emax);
+              else
+                mpfr_nexttozero (a);
               inexact = -MPFR_INT_SIGN (a);
             }
+          else if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+            inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
           MPFR_RET (inexact);
         }
     }
diff --git a/tests/tsub.c b/tests/tsub.c
index a37399d81..e462b780e 100644
--- a/tests/tsub.c
+++ b/tests/tsub.c
@@ -629,6 +629,135 @@ check_rounding (void)
     }
 }
 
+static void
+check_max_almosteven (void)
+{
+  mpfr_exp_t old_emin, old_emax;
+  mpfr_exp_t emin[2] = { MPFR_EMIN_MIN, -1000 };
+  mpfr_exp_t emax[2] = { MPFR_EMAX_MAX, 1000 };
+  int i;
+
+  old_emin = mpfr_get_emin ();
+  old_emax = mpfr_get_emax ();
+
+  for (i = 0; i < 2; i++)
+    {
+      mpfr_t a1, a2, b, c;
+      mpfr_prec_t p;
+      int neg, j, rnd;
+
+      set_emin (emin[i]);
+      set_emax (emax[i]);
+
+      p = MPFR_PREC_MIN + randlimb () % 70;
+      mpfr_init2 (a1, p);
+      mpfr_init2 (a2, p);
+      mpfr_init2 (b, p+1);
+      mpfr_init2 (c, MPFR_PREC_MIN);
+
+      mpfr_setmax (b, 0);
+      mpfr_set_ui (c, 1, MPFR_RNDN);
+
+      for (neg = 0; neg < 2; neg++)
+        {
+          for (j = 1; j >= 0; j--)
+            {
+              mpfr_set_exp (b, __gmpfr_emax - j);
+              RND_LOOP (rnd)
+                {
+                  mpfr_flags_t flags1, flags2;
+                  int inex1, inex2;
+
+                  flags1 = MPFR_FLAGS_INEXACT;
+                  if (rnd == MPFR_RNDN || MPFR_IS_LIKE_RNDZ (rnd, neg))
+                    {
+                      inex1 = neg ? 1 : -1;
+                      mpfr_setmax (a1, __gmpfr_emax - j);
+                    }
+                  else
+                    {
+                      inex1 = neg ? -1 : 1;
+                      if (j == 0)
+                        {
+                          flags1 |= MPFR_FLAGS_OVERFLOW;
+                          mpfr_set_inf (a1, 1);
+                        }
+                      else
+                        {
+                          mpfr_setmin (a1, __gmpfr_emax);
+                        }
+                    }
+                  MPFR_SET_SIGN (a1, neg ? -1 : 1);
+
+                  mpfr_clear_flags ();
+                  inex2 = mpfr_sub (a2, b, c, (mpfr_rnd_t) rnd);
+                  flags2 = __gmpfr_flags;
+
+                  if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                         mpfr_equal_p (a1, a2)))
+                    {
+                      printf ("Error 1 in check_max_almosteven for %s,"
+                              " i = %d, j = %d, neg = %d\n",
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
+                              i, j, neg);
+                      printf ("     b = ");
+                      mpfr_dump (b);
+                      printf ("Expected ");
+                      mpfr_dump (a1);
+                      printf ("  with inex = %d, flags =", inex1);
+                      flags_out (flags1);
+                      printf ("Got      ");
+                      mpfr_dump (a2);
+                      printf ("  with inex = %d, flags =", inex2);
+                      flags_out (flags2);
+                      exit (1);
+                    }
+
+                  if (i == 0)
+                    break;
+
+                  mpfr_clear_flags ();
+                  set_emin (MPFR_EMIN_MIN);
+                  set_emax (MPFR_EMAX_MAX);
+                  inex2 = mpfr_sub (a2, b, c, (mpfr_rnd_t) rnd);
+                  set_emin (emin[i]);
+                  set_emax (emax[i]);
+                  inex2 = mpfr_check_range (a2, inex2, (mpfr_rnd_t) rnd);
+                  flags2 = __gmpfr_flags;
+
+                  if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                         mpfr_equal_p (a1, a2)))
+                    {
+                      printf ("Error 2 in check_max_almosteven for %s,"
+                              " i = %d, j = %d, neg = %d\n",
+                              mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
+                              i, j, neg);
+                      printf ("     b = ");
+                      mpfr_dump (b);
+                      printf ("Expected ");
+                      mpfr_dump (a1);
+                      printf ("  with inex = %d, flags =", inex1);
+                      flags_out (flags1);
+                      printf ("Got      ");
+                      mpfr_dump (a2);
+                      printf ("  with inex = %d, flags =", inex2);
+                      flags_out (flags2);
+                      exit (1);
+                    }
+                }
+            }  /* j */
+
+          mpfr_neg (b, b, MPFR_RNDN);
+          mpfr_neg (c, c, MPFR_RNDN);
+        }  /* neg */
+
+      mpfr_clears (a1, a2, b, c, (mpfr_ptr) 0);
+    }  /* i */
+
+  set_emin (old_emin);
+  set_emax (old_emax);
+}
+
 #define TEST_FUNCTION test_sub
 #define TWO_ARGS
 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
@@ -646,6 +775,7 @@ main (void)
   check_rounding ();
   check_diverse ();
   check_inexact ();
+  check_max_almosteven ();
   bug_ddefour ();
   for (p=2; p<200; p++)
     for (i=0; i<50; i++)