[src/subnormal.c] Set the inexact and underflow flags when needed. Fixes

  https://gforge.inria.fr/tracker/index.php?func=detail&aid=7958&group_id=136&atid=622
[tests/tsubnormal.c] Improved the tests, in particular for the flags.
[doc/mpfr.texi] Documented the change.

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

1 files changed, 39 insertions, 22 deletions

diff --git a/src/subnormal.c b/src/subnormal.c
index 5c13340e2..60e8c428f 100644
--- a/src/subnormal.c
+++ b/src/subnormal.c
@@ -33,33 +33,44 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
      1            |   1         |  ?       | AddOneUlp |
 
    For other rounding mode, there isn't such a problem.
-   Just round it again and merge the inexact flags.
+   Just round it again and merge the ternary values.
+
+   Set the inexact flag if the returned ternary value is non-zero.
+   Set the underflow flag if a second rounding occurred (whether this
+   rounding is exact or not). See
+     http://websympa.loria.fr/wwsympa/arc/mpfr/2009-06/msg00000.html
+     http://websympa.loria.fr/wwsympa/arc/mpfr/2009-06/msg00008.html
+     http://websympa.loria.fr/wwsympa/arc/mpfr/2009-06/msg00010.html
 */
 
 int
 mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
 {
-  int inexact = 0;
+  int sign;
 
-  /* The subnormal exponent range are from:
-      mpfr_emin to mpfr_emin + MPFR_PREC(y) - 1  */
+  /* The subnormal exponent range is [ emin, emin + MPFR_PREC(y) - 2 ] */
   if (MPFR_LIKELY (MPFR_IS_SINGULAR (y)
                    || (MPFR_GET_EXP (y) >=
                        __gmpfr_emin + (mpfr_exp_t) MPFR_PREC (y) - 1)))
-    inexact = old_inexact;
+    MPFR_RET (old_inexact);
+
+  mpfr_set_underflow ();
+  sign = MPFR_SIGN (y);
 
   /* We have to emulate one bit rounding if EXP(y) = emin */
-  else if (MPFR_GET_EXP (y) == __gmpfr_emin)
+  if (MPFR_GET_EXP (y) == __gmpfr_emin)
     {
       /* If this is a power of 2, we don't need rounding.
          It handles cases when rouding away and y=0.1*2^emin */
       if (mpfr_powerof2_raw (y))
-        inexact = old_inexact;
+        MPFR_RET (old_inexact);
+
       /* We keep the same sign for y.
          Assuming Y is the real value and y the approximation
          and since y is not a power of 2:  0.5*2^emin < Y < 1*2^emin
-         We also know the direction of the error thanks to inexact flag */
-      else if (rnd == MPFR_RNDN)
+         We also know the direction of the error thanks to ternary value. */
+
+      if (rnd == MPFR_RNDN)
         {
           mp_limb_t *mant, rb ,sb;
           mp_size_t s;
@@ -77,8 +88,8 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
             goto set_min_p1;
           /* Rounding bit is 1 and sticky bit is 0.
              We need to examine old inexact flag to conclude. */
-          if ((old_inexact > 0 && MPFR_SIGN (y) > 0) ||
-              (old_inexact < 0 && MPFR_SIGN (y) < 0))
+          if ((old_inexact > 0 && sign > 0) ||
+              (old_inexact < 0 && sign < 0))
             goto set_min;
           /* If inexact != 0, return 0.1*2^(emin+1).
              Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0
@@ -91,37 +102,41 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
         {
         set_min:
           mpfr_setmin (y, __gmpfr_emin);
-          inexact = -MPFR_SIGN (y);
+          MPFR_RET (-sign);
         }
       else
         {
         set_min_p1:
           /* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */
           mpfr_setmin (y, __gmpfr_emin + 1);
-          inexact = MPFR_SIGN (y);
+          MPFR_RET (sign);
         }
     }
-
   else /* Hard case: It is more or less the same problem than mpfr_cache */
     {
       mpfr_t dest;
       mpfr_prec_t q;
-      int sign;
+      int inexact, inex2;
+
+      MPFR_ASSERTD (MPFR_GET_EXP (y) > __gmpfr_emin);
 
       /* Compute the intermediary precision */
       q = (mpfr_uexp_t) MPFR_GET_EXP (y) - __gmpfr_emin + 1;
+      MPFR_ASSERTD (q >= MPFR_PREC_MIN && q < MPFR_PREC (y));
+
+      /* TODO: perform the rounding in place. */
       mpfr_init2 (dest, q);
       /* Round y in dest */
-      sign = MPFR_SIGN (y);
       MPFR_SET_EXP (dest, MPFR_GET_EXP (y));
       MPFR_SET_SIGN (dest, sign);
       MPFR_RNDRAW_EVEN (inexact, dest,
                         MPFR_MANT (y), MPFR_PREC (y), rnd, sign,
-                        MPFR_SET_EXP (dest, MPFR_GET_EXP (dest)+1));
+                        MPFR_SET_EXP (dest, MPFR_GET_EXP (dest) + 1));
       if (MPFR_LIKELY (old_inexact != 0))
         {
-          if (MPFR_UNLIKELY(rnd == MPFR_RNDN && (inexact == MPFR_EVEN_INEX
-                                               || inexact == -MPFR_EVEN_INEX)))
+          if (MPFR_UNLIKELY (rnd == MPFR_RNDN &&
+                             (inexact == MPFR_EVEN_INEX ||
+                              inexact == -MPFR_EVEN_INEX)))
             {
               /* if both roundings are in the same direction, we have to go
                  back in the other direction */
@@ -137,10 +152,12 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
           else if (MPFR_UNLIKELY (inexact == 0))
             inexact = old_inexact;
         }
-      old_inexact = mpfr_set (y, dest, rnd);
-      MPFR_ASSERTN (old_inexact == 0);
+
+      inex2 = mpfr_set (y, dest, rnd);
+      MPFR_ASSERTN (inex2 == 0);
       MPFR_ASSERTN (MPFR_IS_PURE_FP (y));
       mpfr_clear (dest);
+
+      MPFR_RET (inexact);
     }
-  return inexact;
 }