more work on the faithful branch

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

8 files changed, 128 insertions, 21 deletions

diff --git a/doc/mpfr.texi b/doc/mpfr.texi
index 6492967ad..fbcf9f17f 100644
--- a/doc/mpfr.texi
+++ b/doc/mpfr.texi
@@ -741,13 +741,14 @@ input and output.)
 @comment  node-name,  next,  previous,  up
 @section Rounding Modes
 
-The following five rounding modes are supported:
+The following rounding modes are supported:
 @itemize @bullet
 @item @code{MPFR_RNDN}: round to nearest (roundTiesToEven in IEEE 754-2008),
 @item @code{MPFR_RNDZ}: round toward zero (roundTowardZero in IEEE 754-2008),
 @item @code{MPFR_RNDU}: round toward plus infinity (roundTowardPositive in IEEE 754-2008),
 @item @code{MPFR_RNDD}: round toward minus infinity (roundTowardNegative in IEEE 754-2008),
 @item @code{MPFR_RNDA}: round away from zero.
+@item @code{MPFR_RNDF}: faithful rounding.
 @end itemize
 
 The @samp{round to nearest} mode works as in the IEEE 754 standard: in
@@ -758,6 +759,15 @@ rounded to (10.0)=2 with a precision of two bits, and not to (11.0)=3.
 This rule avoids the @dfn{drift} phenomenon mentioned by Knuth in volume 2
 of The Art of Computer Programming (Section 4.2.2).
 
+The @code{MPFR_RNDF} mode works as follows: the computed value is either
+that corresponding to @code{MPFR_RNDD} or that corresponding to
+@code{MPFR_RNDU}, except in the case where those values are identical,
+in which case the two adjacent floating-point numbers can also be returned.
+In summary the error is always at most one unit-in-last-place of the computed
+result, which can take two or three possible values.
+For @code{MPFR_RNDF} the ternary value makes no sense.
+Should the flags be correctly set?
+
 @anchor{ternary value}@cindex Ternary value
 Most MPFR functions take as first argument the destination variable, as
 second and following arguments the input variables, as last argument a
diff --git a/tests/tdiv.c b/tests/tdiv.c
index 772ce745a..eb46a004d 100644
--- a/tests/tdiv.c
+++ b/tests/tdiv.c
@@ -29,6 +29,8 @@ check_equal (mpfr_srcptr a, mpfr_srcptr a2, char *s,
   if ((MPFR_IS_NAN (a) && MPFR_IS_NAN (a2)) ||
       mpfr_equal_p (a, a2))
     return;
+  if (r == MPFR_RNDF) /* RNDF might return different values */
+    return;
   printf ("Error in %s\n", mpfr_print_rnd_mode (r));
   printf ("b  = ");
   mpfr_dump (b);
@@ -88,7 +90,7 @@ mpfr_all_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r)
         {
           __gmpfr_flags = oldflags;
           inex2 = mpfr_div_ui (a2, b, mpfr_get_ui (c, MPFR_RNDN), r);
-          MPFR_ASSERTN (SAME_SIGN (inex2, inex));
+          MPFR_ASSERTN (r == MPFR_RNDF || SAME_SIGN (inex2, inex));
           MPFR_ASSERTN (__gmpfr_flags == newflags);
           check_equal (a, a2, "mpfr_div_ui", b, c, r);
         }
@@ -96,7 +98,7 @@ mpfr_all_div (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r)
         {
           __gmpfr_flags = oldflags;
           inex2 = mpfr_div_si (a2, b, mpfr_get_si (c, MPFR_RNDN), r);
-          MPFR_ASSERTN (SAME_SIGN (inex2, inex));
+          MPFR_ASSERTN (r == MPFR_RNDF || SAME_SIGN (inex2, inex));
           MPFR_ASSERTN (__gmpfr_flags == newflags);
           check_equal (a, a2, "mpfr_div_si", b, c, r);
         }
@@ -412,6 +414,9 @@ check_hard (void)
                 {
                   RND_LOOP(rnd)
                     {
+                      if (rnd == MPFR_RNDF)
+                        continue; /* inexact is undefined */
+
                       inex = test_div (q, u, v, (mpfr_rnd_t) rnd);
                       inex2 = get_inexact (q, u, v);
                       if (inex_cmp (inex, inex2))
@@ -676,7 +681,8 @@ check_inexact (void)
               mpfr_set_prec (y, py);
               mpfr_set_prec (z, py + pu);
                 {
-                  rnd = RND_RAND ();
+                  /* inexact is undefined for RNDF */
+                  do rnd = RND_RAND (); while (rnd == MPFR_RNDF);
                   inexact = test_div (y, x, u, rnd);
                   if (mpfr_mul (z, y, u, rnd))
                     {
@@ -975,7 +981,8 @@ consistency (void)
       mpfr_prec_t px, py, pz, p;
       int inex1, inex2;
 
-      rnd = RND_RAND ();
+      /* inex is undefined for RNDF */
+      do rnd = RND_RAND (); while (rnd == MPFR_RNDF);
       px = (randlimb () % 256) + 2;
       py = (randlimb () % 128) + 2;
       pz = (randlimb () % 256) + 2;
@@ -1000,7 +1007,8 @@ consistency (void)
       MPFR_ASSERTN (!MPFR_IS_NAN (z2));
       if (inex1 != inex2 || mpfr_cmp (z1, z2) != 0)
         {
-          printf ("Consistency error for i = %d\n", i);
+          printf ("Consistency error for i = %d (rnd = %s)\n", i,
+                  mpfr_print_rnd_mode (rnd));
           exit (1);
         }
     }
diff --git a/tests/tgeneric.c b/tests/tgeneric.c
index 382c8ce55..fa24fc21c 100644
--- a/tests/tgeneric.c
+++ b/tests/tgeneric.c
@@ -142,11 +142,16 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 #define REDUCE_EMIN mpfr_get_emin ()
 #endif
 
+/* return non-zero if:
+   (1) either x and y are not NaN and are equal
+   (2) x = y = NaN */
+#define EQUAL(x,y) ((mpfr_nan_p (x) && mpfr_nan_p (y)) || mpfr_equal_p (x, y))
+
 static void
 test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
 {
   mpfr_prec_t prec, xprec, yprec;
-  mpfr_t x, y, z, t, w;
+  mpfr_t x, y, z, t, w, yd, yu;
 #ifdef NEED_U
   mpfr_t u;
 #endif
@@ -166,7 +171,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
   old_emin = mpfr_get_emin ();
   old_emax = mpfr_get_emax ();
 
-  mpfr_inits2 (MPFR_PREC_MIN, x, y, z, t, w, (mpfr_ptr) 0);
+  mpfr_inits2 (MPFR_PREC_MIN, x, y, yd, yu, z, t, w, (mpfr_ptr) 0);
 #ifdef NEED_U
   mpfr_init2 (u, MPFR_PREC_MIN);
 #endif
@@ -178,6 +183,8 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
       mpfr_set_prec (t, prec);
       yprec = prec + 10;
       mpfr_set_prec (y, yprec);
+      mpfr_set_prec (yd, yprec);
+      mpfr_set_prec (yu, yprec);
       mpfr_set_prec (w, yprec);
 
 #if defined(TWO_ARGS) || defined(DOUBLE_ARG1) || defined(DOUBLE_ARG2)
@@ -308,9 +315,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
           oemin = mpfr_get_emin ();
           oemax = mpfr_get_emax ();
           
-          do
-            rnd = RND_RAND ();
-          while (rnd == MPFR_RNDF); /* inex makes no sense */
+          rnd = RND_RAND ();
           mpfr_clear_flags ();
 #ifdef DEBUG_TGENERIC
           TGENERIC_INFO (TEST_FUNCTION, MPFR_PREC (y));
@@ -348,6 +353,61 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                           (compare != 0) ^ (mpfr_inexflag_p () == 0));
           ctrt++;
 
+          /* If rnd = RNDF, check that we obtain the same result as
+             RNDD or RNDU, except when RNDD and RNDU are exact */
+          if (rnd == MPFR_RNDF)
+            {
+              int ok;
+#if defined(TWO_ARGS)
+              TEST_FUNCTION (yd, x, u, MPFR_RNDD);
+              TEST_FUNCTION (yu, x, u, MPFR_RNDU);
+#elif defined(DOUBLE_ARG1)
+              d = mpfr_get_d (u, MPFR_RNDD);
+              TEST_FUNCTION (yd, d, x, MPFR_RNDD);
+              d = mpfr_get_d (u, MPFR_RNDU);
+              TEST_FUNCTION (yu, d, x, MPFR_RNDU);
+#elif defined(DOUBLE_ARG2)
+              d = mpfr_get_d (u, MPFR_RNDD);
+              TEST_FUNCTION (yd, x, d, MPFR_RNDD);
+              d = mpfr_get_d (u, MPFR_RNDU);
+              TEST_FUNCTION (yu, x, d, MPFR_RNDU);
+#elif defined(ULONG_ARG1) && defined(ONE_ARG)
+              TEST_FUNCTION (yd, i, MPFR_RNDD);
+              TEST_FUNCTION (yu, i, MPFR_RNDU);
+#elif defined(ULONG_ARG1)
+              TEST_FUNCTION (yd, i, x, MPFR_RNDD);
+              TEST_FUNCTION (yu, i, x, MPFR_RNDU);
+#elif defined(ULONG_ARG2)
+              TEST_FUNCTION (yd, x, i, MPFR_RNDD);
+              TEST_FUNCTION (yu, x, i, MPFR_RNDU);
+#else
+              TEST_FUNCTION (yd, x, MPFR_RNDD);
+              TEST_FUNCTION (yu, x, MPFR_RNDU);
+#endif
+              ok = EQUAL (y, yd) || EQUAL (y, yu);
+              if (compare == 0 && ok == 0)
+                {
+                   mpfr_nextbelow (yd);
+                   mpfr_nextabove (yu);
+                   ok = EQUAL (y, yd) || EQUAL (y, yu);
+                }
+              if (ok == 0)
+                {
+                   printf ("For RNDF, result does not match RNDD nor RNDU\n");
+                   printf ("x = "); mpfr_dump (x);
+#ifdef NEED_U
+                   printf ("u = "); mpfr_dump (u);
+#endif
+#if defined(ULONG_ARG1) || defined(ULONG_ARG2)
+                   printf ("i = %lu\n", i);
+#endif
+                   printf ("yd (RNDD) = "); mpfr_dump (yd);
+                   printf ("yu (RNDU) = "); mpfr_dump (yu);
+                   printf ("y  (RNDF) = "); mpfr_dump (y);
+                   exit (1);
+                }
+            }
+
           /* Tests in a reduced exponent range. */
           {
             mpfr_flags_t oldflags = flags;
@@ -403,7 +463,11 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                     flags = __gmpfr_flags;
                     mpfr_set_emax (oemax);
                     ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
-                    if (flags != ex_flags)
+                    /* for RNDF, this test makes no sense, since RNDF
+                       might return either the maximal floating-point
+                       value or infinity, and the flags might differ in
+                       those two cases. */
+                    if (flags != ex_flags && rnd != MPFR_RNDF)
                       {
                         printf ("tgeneric: error for " MAKE_STR(TEST_FUNCTION)
                                 ", reduced exponent range [%"
@@ -453,7 +517,11 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                     flags = __gmpfr_flags;
                     mpfr_set_emin (oemin);
                     ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
-                    if (flags != ex_flags)
+                    /* for RNDF, this test makes no sense, since RNDF
+                       might return either the maximal floating-point
+                       value or infinity, and the flags might differ in
+                       those two cases. */
+                    if (flags != ex_flags && rnd != MPFR_RNDF)
                       {
                         printf ("tgeneric: error for " MAKE_STR(TEST_FUNCTION)
                                 ", reduced exponent range [%"
@@ -518,7 +586,8 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                 flags = __gmpfr_flags;
                 mpfr_set_emin (oemin);
                 mpfr_set_emax (oemax);
-                if (! (SAME_VAL (w, y) &&
+                /* that test makes no sense for RNDF */
+                if (rnd != MPFR_RNDF && ! (SAME_VAL (w, y) &&
                        SAME_SIGN (inexact, compare) &&
                        flags == oldflags))
                   {
@@ -628,7 +697,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
 #endif
               if (mpfr_erangeflag_p ())
                 goto next_n;
-              if (! mpfr_equal_p (t, z))
+              if (! mpfr_equal_p (t, z) && rnd != MPFR_RNDF)
                 {
                   printf ("tgeneric: results differ for "
                           MAKE_STR(TEST_FUNCTION) " on\n  x = ");
@@ -654,7 +723,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
                 compare = compare + compare2;
               else
                 compare = inexact; /* cannot determine sign(t-f(x)) */
-              if (! SAME_SIGN (inexact, compare))
+              if (! SAME_SIGN (inexact, compare) && rnd != MPFR_RNDF)
                 {
                   printf ("Wrong inexact flag for rnd=%s: expected %d, got %d"
                           "\n", mpfr_print_rnd_mode (rnd), compare, inexact);
@@ -708,7 +777,7 @@ test_generic (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int nmax)
             ctrn, ctrt);
 #endif
 
-  mpfr_clears (x, y, z, t, w, (mpfr_ptr) 0);
+  mpfr_clears (x, y, yd, yu, z, t, w, (mpfr_ptr) 0);
 #ifdef NEED_U
   mpfr_clear (u);
 #endif
diff --git a/tests/tmul.c b/tests/tmul.c
index c503ab2ab..583623963 100644
--- a/tests/tmul.c
+++ b/tests/tmul.c
@@ -227,9 +227,19 @@ check_exact (void)
               printf ("unexpected inexact return value\n");
               exit (1);
             }
-          if ((inexact == 0) && mpfr_cmp (c, d))
+          if ((inexact == 0) && mpfr_cmp (c, d) && rnd != MPFR_RNDF)
             {
-              printf ("inexact=0 but results differ\n");
+              printf ("rnd=%s: inexact=0 but results differ\n",
+                      mpfr_print_rnd_mode (rnd));
+              printf ("a=");
+              mpfr_out_str (stdout, 2, 0, a, rnd);
+              printf ("\nb=");
+              mpfr_out_str (stdout, 2, 0, b, rnd);
+              printf ("\nc=");
+              mpfr_out_str (stdout, 2, 0, c, rnd);
+              printf ("\nd=");
+              mpfr_out_str (stdout, 2, 0, d, rnd);
+              printf ("\n");
               exit (1);
             }
           else if (inexact && (mpfr_cmp (c, d) == 0) && rnd != MPFR_RNDF)
diff --git a/tests/tmul_2exp.c b/tests/tmul_2exp.c
index 5aef54f06..e89179ddc 100644
--- a/tests/tmul_2exp.c
+++ b/tests/tmul_2exp.c
@@ -72,6 +72,9 @@ underflow (mpfr_exp_t e)
             RND_LOOP (rnd)
               for (k = 1; k <= 4; k++)
                 {
+                  if (rnd == MPFR_RNDF)
+                    continue; /* the test below makes no sense for RNDF */
+
                   /* The following one is assumed to be correct. */
                   inex1 = mpfr_mul_2si (y, x, e, MPFR_RNDN);
                   MPFR_ASSERTN (inex1 == 0);
@@ -270,6 +273,9 @@ overflow0 (mpfr_exp_t emax)
           int inex1, inex2;
           mpfr_flags_t flags1, flags2;
 
+          if (r == MPFR_RNDF)
+            continue; /* the test below makes no sense for RNDF */
+
           /* Even if there isn't an overflow (rounding ~ toward zero),
              the result is the same as the one of an overflow. */
           inex1 = mpfr_overflow (y1, (mpfr_rnd_t) r, neg ? -1 : 1);
diff --git a/tests/tmul_ui.c b/tests/tmul_ui.c
index 79f14e960..5a72ebba7 100644
--- a/tests/tmul_ui.c
+++ b/tests/tmul_ui.c
@@ -45,6 +45,9 @@ check_inexact (mpfr_prec_t p)
   for (q = MPFR_PREC_MIN; q <= p; q++)
     for (rnd = 0; rnd < MPFR_RND_MAX; rnd++)
       {
+        if (rnd == MPFR_RNDF)
+          continue; /* inexact is undefined */
+
         mpfr_set_prec (y, q);
         inexact = mpfr_mul_ui (y, x, u, (mpfr_rnd_t) rnd);
         cmp = mpfr_cmp (y, z);
diff --git a/tests/tsqr.c b/tests/tsqr.c
index 199fa4692..ec585d9ce 100644
--- a/tests/tsqr.c
+++ b/tests/tsqr.c
@@ -73,7 +73,8 @@ check_random (mpfr_prec_t p)
             inexact2 = mpfr_sqr (z, x, (mpfr_rnd_t) r);
             if (mpfr_cmp (y, z))
               error1 ((mpfr_rnd_t) r,p,x,y,z);
-            if (inexact_sign (inexact1) != inexact_sign (inexact2))
+            if (inexact_sign (inexact1) != inexact_sign (inexact2) &&
+                r != MPFR_RNDF)
               error2 ((mpfr_rnd_t) r,p,x,y,inexact1,inexact2);
           }
     }
diff --git a/tests/tui_div.c b/tests/tui_div.c
index e5fab1779..f3ef8a549 100644
--- a/tests/tui_div.c
+++ b/tests/tui_div.c
@@ -77,7 +77,7 @@ check_inexact (void)
                   exit (1);
                 }
               cmp = mpfr_cmp_ui (z, u);
-              if (! SAME_SIGN (inexact, cmp))
+              if (rnd != MPFR_RNDF && ! SAME_SIGN (inexact, cmp))
                 {
                   printf ("Wrong inexact flag for u=%lu, rnd=%s\n",
                           u, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));