[src/subnormal.c] experimental support for RNDNA

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

1 files changed, 22 insertions, 6 deletions

diff --git a/src/subnormal.c b/src/subnormal.c
index 70b436e37..be4333203 100644
--- a/src/subnormal.c
+++ b/src/subnormal.c
@@ -70,7 +70,7 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
          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 ternary value. */
 
-      if (rnd == MPFR_RNDN)
+      if (rnd == MPFR_RNDN || rnd == MPFR_RNDNA)
         {
           mp_limb_t *mant, rb ,sb;
           mp_size_t s;
@@ -95,7 +95,7 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
              Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0
              So we have 0.1100000000000000000000000*2^emin exactly.
              We return 0.1*2^(emin+1) according to the even-rounding
-             rule on subnormals. */
+             rule on subnormals. Note the same holds for RNDNA. */
           goto set_min_p1;
         }
       else if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y)))
@@ -129,17 +129,20 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
       /* Round y in dest */
       MPFR_SET_EXP (dest, MPFR_GET_EXP (y));
       MPFR_SET_SIGN (dest, sign);
+      mpfr_rnd_t rnd2 = (rnd == MPFR_RNDNA) ? MPFR_RNDN : rnd;
       MPFR_RNDRAW_EVEN (inexact, dest,
-                        MPFR_MANT (y), MPFR_PREC (y), rnd, sign,
+                        MPFR_MANT (y), MPFR_PREC (y), rnd2, sign,
                         MPFR_SET_EXP (dest, MPFR_GET_EXP (dest) + 1));
       if (MPFR_LIKELY (old_inexact != 0))
         {
-          if (MPFR_UNLIKELY (rnd == MPFR_RNDN &&
+          if (MPFR_UNLIKELY (rnd2 == 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 */
+              /* If both roundings are in the same direction,
+                 we have to go back in the other direction.
+                 For MPFR_RNDNA it is the same, since we are not
+                 exactly in the middle case (old_inexact != 0). */
               if (SAME_SIGN (inexact, old_inexact))
                 {
                   if (SAME_SIGN (inexact, MPFR_INT_SIGN (y)))
@@ -155,6 +158,19 @@ mpfr_subnormalize (mpfr_ptr y, int old_inexact, mpfr_rnd_t rnd)
           else if (MPFR_UNLIKELY (inexact == 0))
             inexact = old_inexact;
         }
+      else if (rnd == MPFR_RNDNA &&
+               (inexact == MPFR_EVEN_INEX || inexact == -MPFR_EVEN_INEX))
+        {
+          /* We are in the middle case: since we used RNDN to round, we should
+             round in the opposite direction when inexact has the opposite
+             sign of y. */
+          if (!SAME_SIGN (inexact, MPFR_INT_SIGN (y)))
+            {
+              mpfr_nexttoinf (dest);
+              MPFR_ASSERTD(!MPFR_IS_INF (dest));
+              inexact = -inexact;
+            }
+        }
 
       inex2 = mpfr_set (y, dest, rnd);
       MPFR_ASSERTN (inex2 == 0);