Avoid integer overflow in MPFR_FAST_COMPUTE_IF_SMALL_INPUT.

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

16 files changed, 64 insertions, 49 deletions

diff --git a/BUGS b/BUGS
index c9f05e0dc..bdc5d573b 100644
--- a/BUGS
+++ b/BUGS
@@ -38,8 +38,7 @@ Known bugs:
   small exponent).
 
 * Incorrect behavior (possible infinite loop, e.g. in mpfr_exp2) in some
-  functions on tiny arguments, e.g. +/- 2^(emin-1), due to an integer
-  overflow in MPFR_FAST_COMPUTE_IF_SMALL_INPUT.
+  functions on tiny arguments, e.g. +/- 2^(emin-1).
 
 * The mpfr_fma function behaves incorrectly if the multiplication overflows
   or underflows. The overflow case has been fixed except in some corner
diff --git a/asin.c b/asin.c
index 12a54d1e0..32fd9866a 100644
--- a/asin.c
+++ b/asin.c
@@ -52,7 +52,8 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
     }
 
   /* asin(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (asin,x, -2*MPFR_GET_EXP (x)+2,1,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (asin, x, -2 * MPFR_GET_EXP (x), 2, 1,
+                                    rnd_mode, {});
 
   /* Set x_p=|x| (x is a normal number) */
   mpfr_init2 (xp, MPFR_PREC (x));
diff --git a/asinh.c b/asinh.c
index 726088a6a..8ca1bc27a 100644
--- a/asinh.c
+++ b/asinh.c
@@ -63,7 +63,8 @@ mpfr_asinh (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
     }
 
   /* asinh(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+2,0,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * MPFR_GET_EXP (x), 2, 0,
+                                    rnd_mode, {});
 
   Ny = MPFR_PREC (y);   /* Precision of output variable */
 
diff --git a/atan.c b/atan.c
index b1d7b99ea..755c89332 100644
--- a/atan.c
+++ b/atan.c
@@ -250,7 +250,8 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
   /* atan(x) = x - x^3/3 + x^5/5...
      so the error is < 2^(3*EXP(x)-1)
      so `EXP(x)-(3*EXP(x)-1)` = -2*EXP(x)+1 */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan,x, -2*MPFR_GET_EXP (x)+1,0,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan, x, -2 * MPFR_GET_EXP (x), 1, 0,
+                                    rnd_mode, {});
 
   /* Set x_p=|x| */
   MPFR_TMP_INIT_ABS (xp, x);
diff --git a/atanh.c b/atanh.c
index d5715555c..7850e1007 100644
--- a/atanh.c
+++ b/atanh.c
@@ -72,7 +72,8 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode)
     }
 
   /* atanh(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP (xt)+1,1,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP (xt), 1, 1,
+                                    rnd_mode, {});
 
   MPFR_SAVE_EXPO_MARK (expo);
 
diff --git a/cos.c b/cos.c
index 6633f61e8..c3e436411 100644
--- a/cos.c
+++ b/cos.c
@@ -203,8 +203,8 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
   MPFR_SAVE_EXPO_MARK (expo);
 
   /* cos(x) = 1-x^2/2 + ..., so error < 2^(2*EXP(x)-1) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, 0-2*MPFR_GET_EXP (x)+1,0,
-                                    rnd_mode, inexact = _inexact; goto end);
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, -2 * MPFR_GET_EXP (x), 1,
+                                    0, rnd_mode, inexact = _inexact; goto end);
 
   /* Compute initial precision */
   precy = MPFR_PREC (y);
diff --git a/erfc.c b/erfc.c
index d692b9fee..998ab49cf 100644
--- a/erfc.c
+++ b/erfc.c
@@ -153,8 +153,8 @@ mpfr_erfc (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
   MPFR_SAVE_EXPO_MARK (expo);
 
   /* erfc(x) ~ 1, with error < 2^(EXP(x)+1) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, 0-MPFR_GET_EXP (x)-1,
-                                    MPFR_SIGN(x) < 0,
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, - MPFR_GET_EXP (x) - 1,
+                                    0, MPFR_SIGN(x) < 0,
                                     rnd, inex = _inexact; goto end);
 
   prec = MPFR_PREC (y) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 3;
diff --git a/expm1.c b/expm1.c
index 0cb33257c..6b908bc84 100644
--- a/expm1.c
+++ b/expm1.c
@@ -70,9 +70,9 @@ mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode)
       /* For -1 < x < 0, abs(expm1(x)-x) < x^2/2.
          For 0 < x < 1,  abs(expm1(x)-x) < x^2. */
       if (MPFR_IS_POS (x))
-        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, rnd_mode, {});
+        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {});
       else
-        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, 1 - ex, 0, rnd_mode, {});
+        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, 0, rnd_mode, {});
     }
 
   if (MPFR_IS_NEG (x) && ex > 5)  /* x <= -32 */
@@ -93,7 +93,7 @@ mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode)
                <= 2^max(MPFR_EMIN_MIN,-LONG_MAX,ceil(x/ln(2)+epsilon))
          with epsilon > 0 */
       mpfr_clear (t);
-      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, minus_one, err, 0, rnd_mode,
+      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, minus_one, err, 0, 0, rnd_mode,
                                         { mpfr_clear (minus_one); });
       mpfr_clear (minus_one);
     }
diff --git a/log1p.c b/log1p.c
index 6582fc5b2..1745c6e1b 100644
--- a/log1p.c
+++ b/log1p.c
@@ -70,9 +70,9 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
       /* For x > 0,    abs(log(1+x)-x) < x^2/2.
          For x > -0.5, abs(log(1+x)-x) < x^2. */
       if (MPFR_IS_POS (x))
-        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex - 1, 0, rnd_mode, {});
+        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex - 1, 0, 0, rnd_mode, {});
       else
-        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, rnd_mode, {});
+        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {});
     }
 
   comp = mpfr_cmp_si (x, -1);
diff --git a/mpfr-impl.h b/mpfr-impl.h
index 4bbbed71b..9e6c0d5b1 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -1138,9 +1138,9 @@ typedef struct {
    x must not be a singular value (NAN, INF or ZERO).
 
    y is the destination (a mpfr_t), x the value to set (a mpfr_t),
-   err the error term (a mp_exp_t), dir (an int) is the direction of
-   the commited error (if dir = 0, it rounds towards 0, if dir=1,
-   it rounds away from 0), rnd the rounding mode.
+   err1+err2 with err2 <= 3 the error term (mp_exp_t's), dir (an int) is
+   the direction of the commited error (if dir = 0, it rounds towards 0,
+   if dir=1, it rounds away from 0), rnd the rounding mode.
 
    It returns from the function a ternary value in case of success.
    If you want to free something, you must fill the "extra" field
@@ -1148,20 +1148,28 @@ typedef struct {
 
    The test is less restrictive thant necessary, but the function
    will finish the check itself.
+
+   Note: err1 + err2 is allowed to overflow as mp_exp_t, but it must give
+   its real value as mpfr_uexp_t.
 */
-#define MPFR_FAST_COMPUTE_IF_SMALL_INPUT(y,x,err,dir,rnd,extra)     \
-  do {                                                              \
-   mp_exp_t _err = (err);                                           \
-   if (MPFR_UNLIKELY (_err > 0                                      \
-                      && (mpfr_uexp_t) _err > MPFR_PREC (y) + 1))   \
-    {                                                               \
-      int _inexact = mpfr_round_near_x ((y),(x),(err),(dir),(rnd)); \
-      if (_inexact != 0)                                            \
-        {                                                           \
-         extra;                                                     \
-         return _inexact;                                           \
-        }                                                           \
-    }                                                               \
+#define MPFR_FAST_COMPUTE_IF_SMALL_INPUT(y,x,err1,err2,dir,rnd,extra)   \
+  do {                                                                  \
+    mpfr_ptr _y = (y);                                                  \
+    mp_exp_t _err1 = (err1);                                            \
+    mp_exp_t _err2 = (err2);                                            \
+    if (_err1 > 0)                                                      \
+      {                                                                 \
+        mpfr_uexp_t _err = (mpfr_uexp_t) _err1 + _err2;                 \
+        if (MPFR_UNLIKELY (_err > MPFR_PREC (_y) + 1))                  \
+          {                                                             \
+            int _inexact = mpfr_round_near_x (_y,(x),_err,(dir),(rnd)); \
+            if (_inexact != 0)                                          \
+              {                                                         \
+                extra;                                                  \
+                return _inexact;                                        \
+              }                                                         \
+          }                                                             \
+      }                                                                 \
   } while (0)
 
 /******************************************************
@@ -1508,7 +1516,8 @@ __MPFR_DECLSPEC void mpfr_dump_mant _MPFR_PROTO ((const mp_limb_t *,
                                                   mp_prec_t));
 
 __MPFR_DECLSPEC int mpfr_round_near_x _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
-                                                    mp_exp_t, int, mp_rnd_t));
+                                                    mpfr_uexp_t, int,
+                                                    mp_rnd_t));
 __MPFR_DECLSPEC void mpfr_abort_prec_max _MPFR_PROTO ((void))
        MPFR_NORETURN_ATTR;
 
diff --git a/round_near_x.c b/round_near_x.c
index 636ecb2e5..c64e60dc7 100644
--- a/round_near_x.c
+++ b/round_near_x.c
@@ -24,7 +24,7 @@ MA 02110-1301, USA. */
 
 /* Use MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */
 
-/* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mp_exp_t err, int dir,
+/* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mpfr_uexp_t err, int dir,
                           mp_rnd_t rnd)
 
    Assuming y = o(f(x)) = o(x + g(x)) with |g(x)| < 2^(EXP(x)-error)
@@ -34,7 +34,7 @@ MA 02110-1301, USA. */
    x must not be a singular value (NAN, INF or ZERO).
 
    y is the destination (a mpfr_t), x the value to set (a mpfr_t),
-   err the error term (a mp_exp_t) such that |g(x)| < 2^(EXP(x)-err),
+   err the error term (a mpfr_uexp_t) such that |g(x)| < 2^(EXP(x)-err),
    dir (an int) is the direction of the error (if dir = 0,
    it rounds towards 0, if dir=1, it rounds away from 0),
    rnd the rounding mode.
@@ -151,7 +151,7 @@ MA 02110-1301, USA. */
  */
 
 int
-mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mp_exp_t err, int dir,
+mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mpfr_uexp_t err, int dir,
                    mp_rnd_t rnd)
 {
   int inexact, sign;
@@ -161,11 +161,14 @@ mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr x, mp_exp_t err, int dir,
   MPFR_ASSERTD (dir == 0 || dir == 1);
 
   /* First check if we can round. The test is more restrictive than
-     necessary. */
-  if (!(err > 0 && (mpfr_uexp_t) err > MPFR_PREC (y) + 1
-        && ((mpfr_uexp_t) err > MPFR_PREC (x)
+     necessary. Note that if err is not representable in an mp_exp_t,
+     then err > MPFR_PREC (x) and the conversion to mp_exp_t will not
+     occur. */
+  if (!(err > MPFR_PREC (y) + 1
+        && (err > MPFR_PREC (x)
             || mpfr_round_p (MPFR_MANT (x), MPFR_LIMB_SIZE (x),
-                             err, MPFR_PREC (y) + (rnd==GMP_RNDN)))))
+                             (mp_exp_t) err,
+                             MPFR_PREC (y) + (rnd == GMP_RNDN)))))
     /* If we assume we can not round, return 0 */
     return 0;
 
diff --git a/sin.c b/sin.c
index e44d215cb..90aa13c36 100644
--- a/sin.c
+++ b/sin.c
@@ -130,12 +130,8 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
     }
 
   /* sin(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
-  /* FIXME: we replaced -2*MPFR_GET_EXP (x)+2 by -2*MPFR_GET_EXP (x)+1,
-     which is still a valid bound, but looses one bit, to overcome the
-     problem when MPFR_GET_EXP(x)=-2^(w-2)+1 on a w-bit machine: then
-     -2*MPFR_GET_EXP (x)+2 = 2^(w-1) overflows. When this problem is
-     fixed in MPFR_FAST_COMPUTE_IF_SMALL_INPUT, replace 1 by 2. */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+1,0,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * MPFR_GET_EXP (x), 2, 0,
+                                    rnd_mode, {});
 
   /* Compute initial precision */
   precy = MPFR_PREC (y);
diff --git a/sinh.c b/sinh.c
index 99ff2a098..b6737ad93 100644
--- a/sinh.c
+++ b/sinh.c
@@ -58,7 +58,8 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
     }
 
   /* sinh(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP(xt)+2,1,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP(xt), 2, 1,
+                                    rnd_mode, {});
 
   MPFR_TMP_INIT_ABS (x, xt);
 
diff --git a/tan.c b/tan.c
index a8e4a3394..ab359784c 100644
--- a/tan.c
+++ b/tan.c
@@ -54,7 +54,8 @@ mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
     }
 
   /* tan(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+1,1,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * MPFR_GET_EXP (x), 1, 1,
+                                    rnd_mode, {});
 
   MPFR_SAVE_EXPO_MARK (expo);
 
diff --git a/tanh.c b/tanh.c
index 4c7dd8313..0a16aae44 100644
--- a/tanh.c
+++ b/tanh.c
@@ -57,7 +57,8 @@ mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode)
     }
 
   /* tanh(x) = x - x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
-  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2*MPFR_GET_EXP(xt)+1,0,rnd_mode,{});
+  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP (xt), 1, 0,
+                                    rnd_mode, {});
 
   MPFR_TMP_INIT_ABS (x, xt);
 
diff --git a/zeta.c b/zeta.c
index caf5ba22a..6bd0c519d 100644
--- a/zeta.c
+++ b/zeta.c
@@ -175,7 +175,8 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
       else
         err = ((mp_exp_t)1) << err;
       err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
-      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 1, rnd_mode,{});
+      MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 0, 1,
+                                        rnd_mode, {});
     }
 
   d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;