[src/set_{si,ui}_2exp.c] Added early underflow/overflow checking to

avoid integer overflow or errors due to special exponent values.

[src/set_uj.c] Added early overflow detection to avoid a possible
integer overflow. Fixed underflow detection.

[tests/tset_si.c] Added tests of mpfr_set_si_2exp and mpfr_set_ui_2exp
in precision 3 with integers from -31 to 31 and exponents MPFR_EXP_MIN,
MPFR_EMIN_MIN-7 to MPFR_EMIN_MIN, MPFR_EMAX_MAX-7 to MPFR_EMAX_MAX,
MPFR_EXP_MAX-7 to MPFR_EXP_MAX (testcases for the above bugs in
src/set_{si,ui}_2exp.c).

[tests/tset_sj.c] Added tests of mpfr_set_sj_2exp and mpfr_set_uj_2exp
in precision 3 with integers from -31 to 31 and the same exponents as
in tset_si.c, but also exponents INTMAX_MIN and INTMAX_MAX (testcases
for the above bugs in src/set_uj.c).

(merged changesets r13808-13809,13826,13832-13835 from the trunk)

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

5 files changed, 311 insertions, 6 deletions

diff --git a/src/set_si_2exp.c b/src/set_si_2exp.c
index e15c9b60f..6b8fa9cd1 100644
--- a/src/set_si_2exp.c
+++ b/src/set_si_2exp.c
@@ -40,6 +40,15 @@ mpfr_set_si_2exp (mpfr_ptr x, long i, mpfr_exp_t e, mpfr_rnd_t rnd_mode)
       mp_limb_t ai, *xp;
       int inex = 0;
 
+      /* Early underflow/overflow checking is necessary to avoid
+         integer overflow or errors due to special exponent values. */
+      if (MPFR_UNLIKELY (e < __gmpfr_emin - (mpfr_exp_t)
+                         (sizeof (unsigned long) * CHAR_BIT + 1)))
+        return mpfr_underflow (x, rnd_mode == MPFR_RNDN ?
+                               MPFR_RNDZ : rnd_mode, i < 0 ? -1 : 1);
+      if (MPFR_UNLIKELY (e >= __gmpfr_emax))
+        return mpfr_overflow (x, rnd_mode, i < 0 ? -1 : 1);
+
       /* FIXME: support int limbs (e.g. 16-bit limbs on 16-bit proc) */
       ai = SAFE_ABS (unsigned long, i);
       MPFR_ASSERTN (SAFE_ABS (unsigned long, i) == ai);
diff --git a/src/set_ui_2exp.c b/src/set_ui_2exp.c
index 172372708..89fc50383 100644
--- a/src/set_ui_2exp.c
+++ b/src/set_ui_2exp.c
@@ -41,6 +41,15 @@ mpfr_set_ui_2exp (mpfr_ptr x, unsigned long i, mpfr_exp_t e, mpfr_rnd_t rnd_mode
       mp_limb_t *xp;
       int inex = 0;
 
+      /* Early underflow/overflow checking is necessary to avoid
+         integer overflow or errors due to special exponent values. */
+      if (MPFR_UNLIKELY (e < __gmpfr_emin - (mpfr_exp_t)
+                         (sizeof (unsigned long) * CHAR_BIT + 1)))
+        return mpfr_underflow (x, rnd_mode == MPFR_RNDN ?
+                               MPFR_RNDZ : rnd_mode, i < 0 ? -1 : 1);
+      if (MPFR_UNLIKELY (e >= __gmpfr_emax))
+        return mpfr_overflow (x, rnd_mode, i < 0 ? -1 : 1);
+
       /* FIXME: support int limbs (e.g. 16-bit limbs on 16-bit proc) */
       MPFR_ASSERTD (i == (mp_limb_t) i);
 
diff --git a/src/set_uj.c b/src/set_uj.c
index 7b94a2a09..5a79ee772 100644
--- a/src/set_uj.c
+++ b/src/set_uj.c
@@ -41,7 +41,7 @@ mpfr_set_uj (mpfr_t x, uintmax_t j, mpfr_rnd_t rnd)
 int
 mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd)
 {
-  int cnt;
+  int cnt, inex;
   mp_size_t i, k;
   mp_limb_t limb;
   mp_limb_t yp[uintmaxpml];
@@ -57,6 +57,10 @@ mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd)
       MPFR_RET(0);
     }
 
+  /* early overflow detection to avoid a possible integer overflow below */
+  if (MPFR_UNLIKELY(e >= __gmpfr_emax))
+    return mpfr_overflow (x, rnd, MPFR_SIGN_POS);
+
   MPFR_ASSERTN (sizeof(uintmax_t) % sizeof(mp_limb_t) == 0);
 
   /* Create an auxiliary var */
@@ -107,7 +111,9 @@ mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd)
   e += k * GMP_NUMB_BITS - cnt;    /* Update Expo */
   MPFR_ASSERTD (MPFR_LIMB_MSB(yp[numberof (yp) - 1]) != 0);
 
-  /* Check expo underflow / overflow (can't use mpfr_check_range) */
+  MPFR_RNDRAW (inex, x, yp, uintmax_bit_size, rnd, MPFR_SIGN_POS, e++);
+
+  /* Check expo underflow / overflow */
   if (MPFR_UNLIKELY(e < __gmpfr_emin))
     {
       /* The following test is necessary because in the rounding to the
@@ -116,16 +122,18 @@ mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd)
        *   _ |x| < 2^(emin-2), or
        *   _ |x| = 2^(emin-2) and the absolute value of the exact
        *     result is <= 2^(emin-2). */
-      if (rnd == MPFR_RNDN && (e+1 < __gmpfr_emin || mpfr_powerof2_raw(y)))
+      if (rnd == MPFR_RNDN &&
+          (e + 1 < __gmpfr_emin ||
+           (mpfr_powerof2_raw (x) && inex >= 0)))
         rnd = MPFR_RNDZ;
       return mpfr_underflow (x, rnd, MPFR_SIGN_POS);
     }
   if (MPFR_UNLIKELY(e > __gmpfr_emax))
     return mpfr_overflow (x, rnd, MPFR_SIGN_POS);
-  MPFR_SET_EXP (y, e);
 
-  /* Final: set x to y (rounding if necessary) */
-  return mpfr_set (x, y, rnd);
+  MPFR_SET_SIGN (x, MPFR_SIGN_POS);
+  MPFR_SET_EXP (x, e);
+  MPFR_RET (inex);
 }
 
 #endif
diff --git a/tests/tset_si.c b/tests/tset_si.c
index fab36efdb..bb194e440 100644
--- a/tests/tset_si.c
+++ b/tests/tset_si.c
@@ -90,6 +90,135 @@ test_2exp (void)
   mpfr_clear (x);
 }
 
+#define REXP 1024
+
+static void
+test_2exp_extreme_aux (void)
+{
+  mpfr_t x1, x2, y;
+  mpfr_exp_t e, ep[1 + 8 * 5], eb[] =
+    { MPFR_EMIN_MIN, -REXP, REXP, MPFR_EMAX_MAX, MPFR_EXP_MAX };
+  mpfr_flags_t flags1, flags2;
+  int i, j, rnd, inex1, inex2;
+  char s;
+
+  ep[0] = MPFR_EXP_MIN;
+  for (i = 0; i < numberof(eb); i++)
+    for (j = 0; j < 8; j++)
+      ep[1 + 8 * i + j] = eb[i] - j;
+
+  mpfr_inits2 (3, x1, x2, (mpfr_ptr) 0);
+  mpfr_init2 (y, 32);
+
+  for (i = 0; i < numberof(ep); i++)
+    for (j = -31; j <= 31; j++)
+      RND_LOOP_NO_RNDF (rnd)
+        {
+          int sign = j < 0 ? -1 : 1;
+
+          /* Compute the expected value, inex and flags */
+          inex1 = mpfr_set_si (y, j, MPFR_RNDN);
+          MPFR_ASSERTN (inex1 == 0);
+          inex1 = mpfr_set (x1, y, (mpfr_rnd_t) rnd);
+          /* x1 is the rounded value and inex1 the ternary value,
+             assuming that the exponent argument is 0 (this is the
+             rounded significand of the final result, assuming an
+             unbounded exponent range). The multiplication by a
+             power of 2 is exact, unless underflow/overflow occurs.
+             The tests on the exponent below avoid integer overflows
+             (ep[i] may take extreme values). */
+          e = mpfr_get_exp (x1);
+          mpfr_clear_flags ();
+          if (j != 0 && ep[i] < __gmpfr_emin - e)  /* underflow */
+            {
+              mpfr_rnd_t r =
+                (rnd == MPFR_RNDN &&
+                 (ep[i] < __gmpfr_emin - mpfr_get_exp (y) - 1 ||
+                  IS_POW2 (sign * j))) ?
+                MPFR_RNDZ : (mpfr_rnd_t) rnd;
+              inex1 = mpfr_underflow (x1, r, sign);
+              flags1 = __gmpfr_flags;
+            }
+          else if (j != 0 && ep[i] > __gmpfr_emax - e)  /* overflow */
+            {
+              inex1 = mpfr_overflow (x1, (mpfr_rnd_t) rnd, sign);
+              flags1 = __gmpfr_flags;
+            }
+          else
+            {
+              if (j != 0)
+                mpfr_set_exp (x1, ep[i] + e);
+              flags1 = inex1 != 0 ? MPFR_FLAGS_INEXACT : 0;
+            }
+
+          /* Test mpfr_set_si_2exp */
+          mpfr_clear_flags ();
+          inex2 = mpfr_set_si_2exp (x2, j, ep[i], (mpfr_rnd_t) rnd);
+          flags2 = __gmpfr_flags;
+
+          if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                 mpfr_equal_p (x1, x2)))
+            {
+              s = 's';
+              goto err_extreme;
+            }
+
+          if (j < 0)
+            continue;
+
+          /* Test mpfr_set_ui_2exp */
+          mpfr_clear_flags ();
+          inex2 = mpfr_set_ui_2exp (x2, j, ep[i], (mpfr_rnd_t) rnd);
+          flags2 = __gmpfr_flags;
+
+          if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                 mpfr_equal_p (x1, x2)))
+            {
+              s = 'u';
+            err_extreme:
+              printf ("Error in extreme mpfr_set_%ci_2exp for i=%d j=%d %s\n",
+                      s, i, j, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+              printf ("emin=%" MPFR_EXP_FSPEC "d "
+                      "emax=%" MPFR_EXP_FSPEC "d\n",
+                      (mpfr_eexp_t) __gmpfr_emin,
+                      (mpfr_eexp_t) __gmpfr_emax);
+              printf ("ep[%d] = %" MPFR_EXP_FSPEC "d\n",
+                      i, (mpfr_eexp_t) ep[i]);
+              printf ("Expected ");
+              mpfr_dump (x1);
+              printf ("with inex = %d and flags =", inex1);
+              flags_out (flags1);
+              printf ("Got      ");
+              mpfr_dump (x2);
+              printf ("with inex = %d and flags =", inex2);
+              flags_out (flags2);
+              exit (1);
+            }
+        }
+
+  mpfr_clears (x1, x2, y, (mpfr_ptr) 0);
+}
+
+static void
+test_2exp_extreme (void)
+{
+  mpfr_exp_t emin, emax;
+
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+
+  set_emin (MPFR_EMIN_MIN);
+  set_emax (MPFR_EMAX_MAX);
+  test_2exp_extreme_aux ();
+
+  set_emin (-REXP);
+  set_emax (REXP);
+  test_2exp_extreme_aux ();
+
+  set_emin (emin);
+  set_emax (emax);
+}
+
 static void
 test_macros (void)
 {
@@ -639,6 +768,7 @@ main (int argc, char *argv[])
   mpfr_clear (x);
 
   test_2exp ();
+  test_2exp_extreme ();
   test_macros ();
   test_macros_keyword ();
   test_get_ui_smallneg ();
diff --git a/tests/tset_sj.c b/tests/tset_sj.c
index a552062b4..e8f81dffc 100644
--- a/tests/tset_sj.c
+++ b/tests/tset_sj.c
@@ -176,6 +176,154 @@ check_set_sj_2exp (void)
   mpfr_clear (x);
 }
 
+#define REXP 1024
+
+static void
+test_2exp_extreme_aux (void)
+{
+  mpfr_t x1, x2, y;
+  mpfr_exp_t e, ep[1 + 8 * 5], eb[] =
+    { MPFR_EMIN_MIN, -REXP, REXP, MPFR_EMAX_MAX, MPFR_EXP_MAX };
+  mpfr_flags_t flags1, flags2;
+  int i, j, rnd, inex1, inex2;
+  char s;
+
+  ep[0] = MPFR_EXP_MIN;
+  for (i = 0; i < numberof(eb); i++)
+    for (j = 0; j < 8; j++)
+      ep[1 + 8 * i + j] = eb[i] - j;
+
+  mpfr_inits2 (3, x1, x2, (mpfr_ptr) 0);
+  mpfr_init2 (y, 32);
+
+  /* The cast to int is needed if numberof(ep) is of unsigned type, in
+     which case the condition without the cast would be always false. */
+  for (i = -2; i < (int) numberof(ep); i++)
+    for (j = -31; j <= 31; j++)
+      RND_LOOP_NO_RNDF (rnd)
+        {
+          int sign = j < 0 ? -1 : 1;
+          intmax_t em;
+
+          if (i < 0)
+            {
+              em = i == -2 ? INTMAX_MIN : INTMAX_MAX;
+              mpfr_clear_flags ();
+              inex1 = j == 0 ?
+                (mpfr_set_zero (x1, +1), 0) : i == -2 ?
+                mpfr_underflow (x1, rnd == MPFR_RNDN ?
+                                MPFR_RNDZ : (mpfr_rnd_t) rnd, sign) :
+                mpfr_overflow (x1, (mpfr_rnd_t) rnd, sign);
+              flags1 = __gmpfr_flags;
+            }
+          else
+            {
+              em = ep[i];
+              /* Compute the expected value, inex and flags */
+              inex1 = mpfr_set_si (y, j, MPFR_RNDN);
+              MPFR_ASSERTN (inex1 == 0);
+              inex1 = mpfr_set (x1, y, (mpfr_rnd_t) rnd);
+              /* x1 is the rounded value and inex1 the ternary value,
+                 assuming that the exponent argument is 0 (this is the
+                 rounded significand of the final result, assuming an
+                 unbounded exponent range). The multiplication by a
+                 power of 2 is exact, unless underflow/overflow occurs.
+                 The tests on the exponent below avoid integer overflows
+                 (ep[i] may take extreme values). */
+              e = mpfr_get_exp (x1);
+              mpfr_clear_flags ();
+              if (j != 0 && ep[i] < __gmpfr_emin - e)  /* underflow */
+                {
+                  mpfr_rnd_t r =
+                    (rnd == MPFR_RNDN &&
+                     (ep[i] < __gmpfr_emin - mpfr_get_exp (y) - 1 ||
+                      IS_POW2 (sign * j))) ?
+                    MPFR_RNDZ : (mpfr_rnd_t) rnd;
+                  inex1 = mpfr_underflow (x1, r, sign);
+                  flags1 = __gmpfr_flags;
+                }
+              else if (j != 0 && ep[i] > __gmpfr_emax - e)  /* overflow */
+                {
+                  inex1 = mpfr_overflow (x1, (mpfr_rnd_t) rnd, sign);
+                  flags1 = __gmpfr_flags;
+                }
+              else
+                {
+                  if (j != 0)
+                    mpfr_set_exp (x1, ep[i] + e);
+                  flags1 = inex1 != 0 ? MPFR_FLAGS_INEXACT : 0;
+                }
+            }
+
+          /* Test mpfr_set_sj_2exp */
+          mpfr_clear_flags ();
+          inex2 = mpfr_set_sj_2exp (x2, j, em, (mpfr_rnd_t) rnd);
+          flags2 = __gmpfr_flags;
+
+          if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                 mpfr_equal_p (x1, x2)))
+            {
+              s = 's';
+              goto err_extreme;
+            }
+
+          if (j < 0)
+            continue;
+
+          /* Test mpfr_set_uj_2exp */
+          mpfr_clear_flags ();
+          inex2 = mpfr_set_uj_2exp (x2, j, em, (mpfr_rnd_t) rnd);
+          flags2 = __gmpfr_flags;
+
+          if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+                 mpfr_equal_p (x1, x2)))
+            {
+              s = 'u';
+            err_extreme:
+              printf ("Error in extreme mpfr_set_%cj_2exp for i=%d j=%d %s\n",
+                      s, i, j, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+              printf ("emin=%" MPFR_EXP_FSPEC "d "
+                      "emax=%" MPFR_EXP_FSPEC "d\n",
+                      (mpfr_eexp_t) __gmpfr_emin,
+                      (mpfr_eexp_t) __gmpfr_emax);
+#ifndef NPRINTF_J
+              printf ("e = %jd\n", em);
+#endif
+              printf ("Expected ");
+              mpfr_dump (x1);
+              printf ("with inex = %d and flags =", inex1);
+              flags_out (flags1);
+              printf ("Got      ");
+              mpfr_dump (x2);
+              printf ("with inex = %d and flags =", inex2);
+              flags_out (flags2);
+              exit (1);
+            }
+        }
+
+  mpfr_clears (x1, x2, y, (mpfr_ptr) 0);
+}
+
+static void
+test_2exp_extreme (void)
+{
+  mpfr_exp_t emin, emax;
+
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+
+  set_emin (MPFR_EMIN_MIN);
+  set_emax (MPFR_EMAX_MAX);
+  test_2exp_extreme_aux ();
+
+  set_emin (-REXP);
+  set_emax (REXP);
+  test_2exp_extreme_aux ();
+
+  set_emin (emin);
+  set_emax (emax);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -185,6 +333,7 @@ main (int argc, char *argv[])
   check_set_uj_2exp ();
   check_set_sj ();
   check_set_sj_2exp ();
+  test_2exp_extreme ();
 
   tests_end_mpfr ();
   return 0;