Fixed mpfr_div_2si, mpfr_div_2ui and mpfr_mul_2si in some corner cases.

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

4 files changed, 54 insertions, 46 deletions

diff --git a/src/div_2si.c b/src/div_2si.c
index 99ef35365..f55057390 100644
--- a/src/div_2si.c
+++ b/src/div_2si.c
@@ -30,11 +30,12 @@ mpfr_div_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd_mode)
   MPFR_LOG_FUNC (("x[%#R]=%R n=%ld rnd=%d", x, x, n, rnd_mode),
                  ("y[%#R]=%R inexact=%d", y, y, inexact));
 
-  inexact = MPFR_UNLIKELY(y != x) ? mpfr_set (y, x, rnd_mode) : 0;
-
-  if (MPFR_LIKELY( MPFR_IS_PURE_FP(y) ))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+    return mpfr_set (y, x, rnd_mode);
+  else
     {
-      mpfr_exp_t exp = MPFR_GET_EXP (y);
+      mpfr_exp_t exp = MPFR_GET_EXP (x);
+      MPFR_SETRAW (inexact, y, x, exp, rnd_mode);
       if (MPFR_UNLIKELY( n > 0 && (__gmpfr_emin > MPFR_EMAX_MAX - n ||
                                    exp < __gmpfr_emin + n)) )
         {
@@ -45,13 +46,12 @@ mpfr_div_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd_mode)
             rnd_mode = MPFR_RNDZ;
           return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
         }
-
-      if (MPFR_UNLIKELY(n < 0 && (__gmpfr_emax < MPFR_EMIN_MIN - n ||
-                                  exp > __gmpfr_emax + n)) )
+      else if (MPFR_UNLIKELY(n < 0 && (__gmpfr_emax < MPFR_EMIN_MIN - n ||
+                                       exp > __gmpfr_emax + n)) )
         return mpfr_overflow (y, rnd_mode, MPFR_SIGN(y));
 
       MPFR_SET_EXP (y, exp - n);
     }
 
-  return inexact;
+  MPFR_RET (inexact);
 }
diff --git a/src/div_2ui.c b/src/div_2ui.c
index 662da208c..4621470a1 100644
--- a/src/div_2ui.c
+++ b/src/div_2ui.c
@@ -30,42 +30,30 @@ mpfr_div_2ui (mpfr_ptr y, mpfr_srcptr x, unsigned long n, mpfr_rnd_t rnd_mode)
   MPFR_LOG_FUNC (("x[%#R]=%R n=%lu rnd=%d", x, x, n, rnd_mode),
                  ("y[%#R]=%R inexact=%d", y, y, inexact));
 
-  /* Most of the times, this function is called with y==x */
-  inexact = MPFR_UNLIKELY(y != x) ? mpfr_set (y, x, rnd_mode) : 0;
-
-  if (MPFR_LIKELY( MPFR_IS_PURE_FP(y)) )
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+    return mpfr_set (y, x, rnd_mode);
+  else
     {
-      /* n will have to be casted to long to make sure that the addition
-         and subtraction below (for overflow detection) are signed */
-      while (MPFR_UNLIKELY(n > LONG_MAX))
-        {
-          int inex2;
+      mpfr_exp_t exp = MPFR_GET_EXP (x);
+      mpfr_uexp_t diffexp;
 
+      MPFR_SETRAW (inexact, y, x, exp, rnd_mode);
+      diffexp = (mpfr_uexp_t) exp - (mpfr_uexp_t) (__gmpfr_emin - 1);
+      if (MPFR_UNLIKELY (n >= diffexp))  /* exp - n <= emin - 1 */
+        {
+          if (rnd_mode == MPFR_RNDN &&
+              (n > diffexp || (inexact >= 0 && mpfr_powerof2_raw (y))))
+            rnd_mode = MPFR_RNDZ;
+          return mpfr_underflow (y, rnd_mode, MPFR_SIGN (y));
+        }
+      /* exp - n >= emin (no underflow, no integer overflow) */
+      while (n > LONG_MAX)
+        {
           n -= LONG_MAX;
-          inex2 = mpfr_div_2ui(y, y, LONG_MAX, rnd_mode);
-          if (inex2)
-            return inex2; /* underflow */
+          exp -= LONG_MAX;  /* note: signed values */
         }
-
-      /* MPFR_EMAX_MAX - (long) n is signed and doesn't lead to an integer
-         overflow; the first test useful so that the real test can't lead
-         to an integer overflow. */
-      {
-        mpfr_exp_t exp = MPFR_GET_EXP (y);
-        if (MPFR_UNLIKELY( __gmpfr_emin > MPFR_EMAX_MAX - (long) n ||
-                           exp < __gmpfr_emin + (long) n) )
-          {
-            if (rnd_mode == MPFR_RNDN &&
-                (__gmpfr_emin > MPFR_EMAX_MAX - (long) (n - 1) ||
-                 exp < __gmpfr_emin + (long) (n - 1) ||
-                 (inexact >= 0 && mpfr_powerof2_raw (y))))
-              rnd_mode = MPFR_RNDZ;
-            return mpfr_underflow (y, rnd_mode, MPFR_SIGN(y));
-          }
-
-        MPFR_SET_EXP(y, exp - (long) n);
-      }
+      MPFR_SET_EXP (y, exp - (long) n);
     }
 
-  return inexact;
+  MPFR_RET (inexact);
 }
diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h
index 09c411ade..44f7419b2 100644
--- a/src/mpfr-impl.h
+++ b/src/mpfr-impl.h
@@ -1325,6 +1325,26 @@ typedef struct {
  (!MPFR_IS_SINGULAR (b) && mpfr_round_p (MPFR_MANT (b), MPFR_LIMB_SIZE (b),  \
                                          (err), (prec) + ((rnd)==MPFR_RNDN)))
 
+/* Copy the sign and the significand, and handle the exponent in exp. */
+#define MPFR_SETRAW(inexact,dest,src,exp,rnd)                           \
+  if (MPFR_UNLIKELY (dest != src))                                      \
+    {                                                                   \
+      MPFR_SET_SIGN (dest, MPFR_SIGN (src));                            \
+      if (MPFR_LIKELY (MPFR_PREC (dest) == MPFR_PREC (src)))            \
+        {                                                               \
+          MPN_COPY (MPFR_MANT (dest), MPFR_MANT (src),                  \
+                    (MPFR_PREC (src) + GMP_NUMB_BITS-1)/GMP_NUMB_BITS); \
+          inexact = 0;                                                  \
+        }                                                               \
+      else                                                              \
+        {                                                               \
+          MPFR_RNDRAW (inexact, dest, MPFR_MANT (src), MPFR_PREC (src), \
+                       rnd, MPFR_SIGN (src), exp++);                    \
+        }                                                               \
+    }                                                                   \
+  else                                                                  \
+    inexact = 0;
+
 /* TODO: fix this description (see round_near_x.c). */
 /* Assuming that the function has a Taylor expansion which looks like:
     y=o(f(x)) = o(v + g(x)) with |g(x)| <= 2^(EXP(v)-err)
diff --git a/src/mul_2si.c b/src/mul_2si.c
index 0ac0ea35f..79e3fa2d9 100644
--- a/src/mul_2si.c
+++ b/src/mul_2si.c
@@ -30,15 +30,15 @@ mpfr_mul_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd_mode)
   MPFR_LOG_FUNC (("x[%#R]=%R n=%ld rnd=%d", x, x, n, rnd_mode),
                  ("y[%#R]=%R inexact=%d", y, y, inexact));
 
-  inexact = MPFR_UNLIKELY(y != x) ? mpfr_set (y, x, rnd_mode) : 0;
-
-  if (MPFR_LIKELY( MPFR_IS_PURE_FP(y)) )
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+    return mpfr_set (y, x, rnd_mode);
+  else
     {
-      mpfr_exp_t exp = MPFR_GET_EXP (y);
+      mpfr_exp_t exp = MPFR_GET_EXP (x);
+      MPFR_SETRAW (inexact, y, x, exp, rnd_mode);
       if (MPFR_UNLIKELY( n > 0 && (__gmpfr_emax < MPFR_EMIN_MIN + n ||
                                    exp > __gmpfr_emax - n)))
         return mpfr_overflow (y, rnd_mode, MPFR_SIGN(y));
-
       else if (MPFR_UNLIKELY(n < 0 && (__gmpfr_emin > MPFR_EMAX_MAX + n ||
                                        exp < __gmpfr_emin - n)))
         {
@@ -52,5 +52,5 @@ mpfr_mul_2si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd_mode)
       MPFR_SET_EXP (y, exp + n);
     }
 
-  return inexact;
+  MPFR_RET (inexact);
 }