Added mpfr_rint_roundeven and mpfr_roundeven functions, with

documentation and tests.

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

5 files changed, 128 insertions, 42 deletions

diff --git a/NEWS b/NEWS
index eaa9f4808..6d3a8cc94 100644
--- a/NEWS
+++ b/NEWS
@@ -33,6 +33,9 @@ Changes from versions 3.1.* to version 3.2.0:
   mpfr_flags_save and mpfr_flags_restore to operate on groups of flags.
 - New functions mpfr_set_float128 and mpfr_get_float128 to convert from/to
   the __float128 type (requires --enable-float128 and compiler support).
+- New functions mpfr_rint_roundeven and mpfr_roundeven, completing the
+  other similar round-to-integer functions for rounding to nearest with
+  the even-rounding rule.
 - New macro mpfr_round_nearest_away to add partial emulation of the
   rounding to nearest-away (as defined in IEEE 754-2008).
 - New functions mpfr_nrandom and mpfr_erandom to generate random numbers
diff --git a/doc/mpfr.texi b/doc/mpfr.texi
index c8068ddf0..be014b28a 100644
--- a/doc/mpfr.texi
+++ b/doc/mpfr.texi
@@ -2573,15 +2573,26 @@ the null-terminator, or a negative value if an error occurred.
 @deftypefunx int mpfr_ceil (mpfr_t @var{rop}, mpfr_t @var{op})
 @deftypefunx int mpfr_floor (mpfr_t @var{rop}, mpfr_t @var{op})
 @deftypefunx int mpfr_round (mpfr_t @var{rop}, mpfr_t @var{op})
+@deftypefunx int mpfr_roundeven (mpfr_t @var{rop}, mpfr_t @var{op})
 @deftypefunx int mpfr_trunc (mpfr_t @var{rop}, mpfr_t @var{op})
 Set @var{rop} to @var{op} rounded to an integer.
 @code{mpfr_rint} rounds to the nearest representable integer in the
-given direction @var{rnd}, @code{mpfr_ceil} rounds
-to the next higher or equal representable integer, @code{mpfr_floor} to
-the next lower or equal representable integer, @code{mpfr_round} to the
-nearest representable integer, rounding halfway cases away from zero
-(as in the roundTiesToAway mode of IEEE 754-2008),
-and @code{mpfr_trunc} to the next representable integer toward zero.
+given direction @var{rnd}, and the other five functions behave in a
+similar way with some fixed rounding mode:
+@itemize @bullet
+@item @code{mpfr_ceil}: to the next higher or equal representable integer
+      (like @code{mpfr_rint} with @code{MPFR_RNDU});
+@item @code{mpfr_floor} to the next lower or equal representable integer
+      (like @code{mpfr_rint} with @code{MPFR_RNDD});
+@item @code{mpfr_round} to the nearest representable integer,
+      rounding halfway cases away from zero
+      (as in the roundTiesToAway mode of IEEE 754-2008);
+@item @code{mpfr_roundeven} to the nearest representable integer,
+      rounding halfway cases with the even-rounding rule
+      (like @code{mpfr_rint} with @code{MPFR_RNDN});
+@item @code{mpfr_trunc} to the next representable integer toward zero
+      (like @code{mpfr_rint} with @code{MPFR_RNDZ}).
+@end itemize
 When @var{op} is a zero or an infinity, set @var{rop} to the same value
 (with the same sign).
 
@@ -2597,14 +2608,10 @@ the inexact flag is set when @var{rop} differs from @var{op}, following
 the ISO C99 rule for the @code{rint} function. If you want the behavior to
 be more like IEEE 754 / ISO TS 18661-1, i.e., the usual behavior where the
 round-to-integer function is regarded as any other mathematical function,
-you should use one the @code{mpfr_rint_*} functions instead (however it is
-not possible to round to nearest with the even rounding rule yet).
-
-Note that @code{mpfr_round} is different from @code{mpfr_rint} called with
-the rounding to nearest mode (where halfway cases are rounded to an even
-integer or significand). Note also that no double rounding is performed; for
-instance, 10.5 (1010.1 in binary) is rounded by @code{mpfr_rint} with
-rounding to nearest to 12 (1100
+you should use one the @code{mpfr_rint_*} functions instead.
+
+Note that no double rounding is performed; for instance, 10.5 (1010.1 in
+binary) is rounded by @code{mpfr_rint} with rounding to nearest to 12 (1100
 in binary) in 2-bit precision, because the two enclosing numbers representable
 on two bits are 8 and 12, and the closest is 12.
 (If one first rounded to an integer, one would round 10.5 to 10 with
@@ -2614,12 +2621,18 @@ even rounding, and then 10 would be rounded to 8 again with even rounding.)
 @deftypefun int mpfr_rint_ceil (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_rint_floor (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_rint_round (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
+@deftypefunx int mpfr_rint_roundeven (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
 @deftypefunx int mpfr_rint_trunc (mpfr_t @var{rop}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd})
-Set @var{rop} to @var{op} rounded to an integer.
-@code{mpfr_rint_ceil} rounds to the next higher or equal integer,
-@code{mpfr_rint_floor} to the next lower or equal integer,
-@code{mpfr_rint_round} to the nearest integer, rounding halfway cases away
-from zero, and @code{mpfr_rint_trunc} to the next integer toward zero.
+Set @var{rop} to @var{op} rounded to an integer:
+@itemize @bullet
+@item @code{mpfr_rint_ceil}: to the next higher or equal integer;
+@item @code{mpfr_rint_floor}: to the next lower or equal integer;
+@item @code{mpfr_rint_round}: to the nearest integer,
+      rounding halfway cases away from zero;
+@item @code{mpfr_rint_roundeven}: to the nearest integer,
+      rounding halfway cases to the nearest even integer;
+@item @code{mpfr_rint_trunc} to the next integer toward zero.
+@end itemize
 If the result is not representable, it is rounded in the direction @var{rnd}.
 When @var{op} is a zero or an infinity, set @var{rop} to the same value
 (with the same sign).
@@ -3636,6 +3649,8 @@ use @code{mpfr_get_z_exp}.
 
 @item @code{mpfr_remainder} and @code{mpfr_remquo} in MPFR 2.3.
 
+@item @code{mpfr_rint_roundeven} and @code{mpfr_roundeven} in MPFR 3.2.
+
 @item @code{mpfr_round_nearest_away} in MPFR 3.2.
 
 @item @code{mpfr_set_divby0} in MPFR 3.1 (new divide-by-zero exception).
diff --git a/src/mpfr.h b/src/mpfr.h
index 7bb8f1307..b3acb929a 100644
--- a/src/mpfr.h
+++ b/src/mpfr.h
@@ -649,10 +649,13 @@ __MPFR_DECLSPEC int mpfr_div_2si _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
                                                long, mpfr_rnd_t));
 
 __MPFR_DECLSPEC int mpfr_rint _MPFR_PROTO((mpfr_ptr,mpfr_srcptr, mpfr_rnd_t));
+__MPFR_DECLSPEC int mpfr_roundeven _MPFR_PROTO((mpfr_ptr, mpfr_srcptr));
 __MPFR_DECLSPEC int mpfr_round _MPFR_PROTO((mpfr_ptr, mpfr_srcptr));
 __MPFR_DECLSPEC int mpfr_trunc _MPFR_PROTO((mpfr_ptr, mpfr_srcptr));
 __MPFR_DECLSPEC int mpfr_ceil _MPFR_PROTO((mpfr_ptr, mpfr_srcptr));
 __MPFR_DECLSPEC int mpfr_floor _MPFR_PROTO((mpfr_ptr, mpfr_srcptr));
+__MPFR_DECLSPEC int mpfr_rint_roundeven _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
+                                                      mpfr_rnd_t));
 __MPFR_DECLSPEC int mpfr_rint_round _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
                                                   mpfr_rnd_t));
 __MPFR_DECLSPEC int mpfr_rint_trunc _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
diff --git a/src/rint.c b/src/rint.c
index ed0494ff5..5969d8853 100644
--- a/src/rint.c
+++ b/src/rint.c
@@ -304,6 +304,14 @@ mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
     }  /* exp > 0, |u| >= 1 */
 }
 
+#undef mpfr_roundeven
+
+int
+mpfr_roundeven (mpfr_ptr r, mpfr_srcptr u)
+{
+  return mpfr_rint (r, u, MPFR_RNDN);
+}
+
 #undef mpfr_round
 
 int
@@ -341,6 +349,32 @@ mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
  * the inexact flag when the argument is not an integer.
  */
 
+#undef mpfr_rint_roundeven
+
+int
+mpfr_rint_roundeven (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
+{
+  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
+    return mpfr_set (r, u, rnd_mode);
+  else
+    {
+      mpfr_t tmp;
+      int inex;
+      mpfr_flags_t saved_flags = __gmpfr_flags;
+      MPFR_BLOCK_DECL (flags);
+
+      mpfr_init2 (tmp, MPFR_PREC (u));
+      /* round(u) is representable in tmp unless an overflow occurs */
+      MPFR_BLOCK (flags, mpfr_roundeven (tmp, u));
+      __gmpfr_flags = saved_flags;
+      inex = (MPFR_OVERFLOW (flags)
+              ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u))
+              : mpfr_set (r, tmp, rnd_mode));
+      mpfr_clear (tmp);
+      return inex;
+    }
+}
+
 #undef mpfr_rint_round
 
 int
diff --git a/tests/trint.c b/tests/trint.c
index fbfa93164..9c9a645a9 100644
--- a/tests/trint.c
+++ b/tests/trint.c
@@ -269,15 +269,18 @@ basic_tests (void)
       for (s = 1; s >= -1; s -= 2)
         for (i = 1; i <= 72; i++)
           {
-            int k, t, u, v, f, e;
+            int k, t, u, v, f, e, b;
 
             for (t = i/4, k = 0; t >= 1 << prec; t >>= 1, k++)
               ;
+            b = !(t & 1);
             t <<= k;
             for (u = (i+3)/4, k = 0; u >= 1 << prec; u = (u+1)/2, k++)
               ;
             u <<= k;
             v = i < (t+u) << 1 ? t : u;
+            if (b)
+              b = i == (t+u) << 1;
             f = t == u ? 0 : i % 4 == 0 ? 1 : 2;
 
             mpfr_set_si_2exp (x, s * i, -2, MPFR_RNDN);
@@ -288,11 +291,14 @@ basic_tests (void)
                 BASIC_TEST (floor, s > 0 ? i/4 : - ((i+3)/4));
                 BASIC_TEST (ceil, s > 0 ? (i+3)/4 : - (i/4));
                 BASIC_TEST (round, s * ((i+2)/4));
+                BASIC_TEST (roundeven, s * (i % 8 == 2 ? i/4 : (i+2)/4));
               }
             BASIC_TEST2 (trunc, s * t, - s * f);
             BASIC_TEST2 (floor, s > 0 ? t : - u, - f);
             BASIC_TEST2 (ceil, s > 0 ? u : - t, f);
             BASIC_TEST2 (round, s * v, v == t ? - s * f : s * f);
+            BASIC_TEST2 (roundeven, s * (b ? t : v),
+                         b || v == t ? - s * f : s * f);
           }
       mpfr_clears (y, z, (mpfr_ptr) 0);
     }
@@ -463,6 +469,12 @@ coverage_03032011 (void)
 #define test_generic test_generic_round
 #include "tgeneric.c"
 
+#define TEST_FUNCTION mpfr_rint_roundeven
+#define TEST_RANDOM_EMIN -20
+#define TEST_RANDOM_ALWAYS_SCALE 1
+#define test_generic test_generic_roundeven
+#include "tgeneric.c"
+
 int
 main (int argc, char *argv[])
 {
@@ -508,30 +520,43 @@ main (int argc, char *argv[])
           mpfr_set_prec (y, p);
           mpfr_set_prec (v, p);
           for (r = 0; r < MPFR_RND_MAX ; r++)
-            for (trint = 0; trint < 3; trint++)
+            for (trint = 0; trint < 4; trint++)
               {
                 if (trint == 2)
                   inexact = mpfr_rint (y, x, (mpfr_rnd_t) r);
+                else if (trint == 3)
+                  {
+                    if (r != MPFR_RNDN)
+                      continue;
+                    inexact = mpfr_round (y, x);
+                  }
                 else if (r == MPFR_RNDN)
-                  inexact = mpfr_round (y, x);
+                  inexact = (trint ? mpfr_roundeven (y, x) :
+                             mpfr_rint_roundeven (y, x, MPFR_RNDZ));
                 else if (r == MPFR_RNDZ)
                   inexact = (trint ? mpfr_trunc (y, x) :
                              mpfr_rint_trunc (y, x, MPFR_RNDZ));
                 else if (r == MPFR_RNDU)
                   inexact = (trint ? mpfr_ceil (y, x) :
                              mpfr_rint_ceil (y, x, MPFR_RNDU));
-                else /* r = MPFR_RNDD */
+                else if (r == MPFR_RNDD)
                   inexact = (trint ? mpfr_floor (y, x) :
                              mpfr_rint_floor (y, x, MPFR_RNDD));
+                else
+                  {
+                    MPFR_ASSERTN (r == MPFR_RNDA);
+                    continue;
+                  }
                 if (mpfr_sub (t, y, x, MPFR_RNDN))
-                  err ("subtraction 1 should be exact",
-                       s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                  err ("subtraction 1 should be exact", s, x, y, p,
+                       (mpfr_rnd_t) r, trint, inexact);
                 sign_t = mpfr_cmp_ui (t, 0);
                 if (trint != 0 &&
                     (((inexact == 0) && (sign_t != 0)) ||
                      ((inexact < 0) && (sign_t >= 0)) ||
                      ((inexact > 0) && (sign_t <= 0))))
-                  err ("wrong inexact flag", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                  err ("wrong inexact flag", s, x, y, p,
+                       (mpfr_rnd_t) r, trint, inexact);
                 if (inexact == 0)
                   continue; /* end of the test for exact results */
 
@@ -539,28 +564,30 @@ main (int argc, char *argv[])
                      && inexact > 0) ||
                     ((r == MPFR_RNDU || (r == MPFR_RNDZ && MPFR_IS_NEG (x)))
                      && inexact < 0))
-                  err ("wrong rounding direction",
-                       s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                  err ("wrong rounding direction", s, x, y, p,
+                       (mpfr_rnd_t) r, trint, inexact);
                 if (inexact < 0)
                   {
                     mpfr_add_ui (v, y, 1, MPFR_RNDU);
                     if (mpfr_cmp (v, x) <= 0)
-                      err ("representable integer between x and its "
-                           "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                      err ("representable integer between x and "
+                           "its rounded value", s, x, y, p,
+                           (mpfr_rnd_t) r, trint, inexact);
                   }
                 else
                   {
                     mpfr_sub_ui (v, y, 1, MPFR_RNDD);
                     if (mpfr_cmp (v, x) >= 0)
-                      err ("representable integer between x and its "
-                           "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                      err ("representable integer between x and "
+                           "its rounded value", s, x, y, p,
+                           (mpfr_rnd_t) r, trint, inexact);
                   }
-                if (r == MPFR_RNDN)
+                if (r == MPFR_RNDN && trint != 0)
                   {
                     int cmp;
                     if (mpfr_sub (u, v, x, MPFR_RNDN))
-                      err ("subtraction 2 should be exact",
-                           s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                      err ("subtraction 2 should be exact", s, x, y, p,
+                           (mpfr_rnd_t) r, trint, inexact);
                     cmp = mpfr_cmp_abs (t, u);
                     if (cmp > 0)
                       err ("faithful rounding, but not the nearest integer",
@@ -575,8 +602,9 @@ main (int argc, char *argv[])
                            mode: round to an even integer or significand. */
                         mpfr_div_2ui (y, y, 1, MPFR_RNDZ);
                         if (!mpfr_integer_p (y))
-                          err ("halfway case for mpfr_rint, result isn't an"
-                               " even integer", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                          err ("halfway case for mpfr_rint, result isn't "
+                               "an even integer", s, x, y, p,
+                               (mpfr_rnd_t) r, trint, inexact);
                         /* If floor(x) and ceil(x) aren't both representable
                            integers, the significand must be even. */
                         mpfr_sub (v, v, y, MPFR_RNDN);
@@ -586,17 +614,19 @@ main (int argc, char *argv[])
                             mpfr_div_2si (y, y, MPFR_EXP (y) - MPFR_PREC (y)
                                           + 1, MPFR_RNDN);
                             if (!mpfr_integer_p (y))
-                              err ("halfway case for mpfr_rint, significand isn't"
-                                   " even", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                              err ("halfway case for mpfr_rint, "
+                                   "significand isn't even", s, x, y, p,
+                                   (mpfr_rnd_t) r, trint, inexact);
                           }
                       }
-                    else
+                    else if (trint == 3)
                       { /* halfway case for mpfr_round: x must have been
                            rounded away from zero. */
                         if ((MPFR_IS_POS (x) && inexact < 0) ||
                             (MPFR_IS_NEG (x) && inexact > 0))
-                          err ("halfway case for mpfr_round, bad rounding"
-                               " direction", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
+                          err ("halfway case for mpfr_round, "
+                               "bad rounding direction", s, x, y, p,
+                               (mpfr_rnd_t) r, trint, inexact);
                       }
                   }
               }
@@ -617,6 +647,7 @@ main (int argc, char *argv[])
   test_generic_floor (2, 300, 20);
   test_generic_ceil (2, 300, 20);
   test_generic_round (2, 300, 20);
+  test_generic_roundeven (2, 300, 20);
 
 #if __MPFR_STDC (199901L)
   if (argc > 1 && strcmp (argv[1], "-s") == 0)