From e0263aed74458d29a8b976cd1491c3bd8641baa4 Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 20 Jul 2015 15:59:51 +0000 Subject: 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 --- NEWS | 3 +++ doc/mpfr.texi | 53 +++++++++++++++++++++++++--------------- src/mpfr.h | 3 +++ src/rint.c | 34 ++++++++++++++++++++++++++ tests/trint.c | 77 +++++++++++++++++++++++++++++++++++++++++------------------ 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) -- cgit v1.2.1