/* Test file for mpfr_rint, mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round, mpfr_rint_trunc, mpfr_rint_floor, mpfr_rint_ceil, mpfr_rint_round. Copyright 2002-2017 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. The GNU MPFR Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU MPFR Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ #include "mpfr-test.h" #if __MPFR_STDC (199901L) # include #endif static void special (void) { mpfr_t x, y; mpfr_exp_t emax; mpfr_init (x); mpfr_init (y); mpfr_set_nan (x); mpfr_rint (y, x, MPFR_RNDN); MPFR_ASSERTN(mpfr_nan_p (y)); mpfr_set_inf (x, 1); mpfr_rint (y, x, MPFR_RNDN); MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0); mpfr_set_inf (x, -1); mpfr_rint (y, x, MPFR_RNDN); MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) < 0); mpfr_set_ui (x, 0, MPFR_RNDN); mpfr_rint (y, x, MPFR_RNDN); MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y)); mpfr_set_ui (x, 0, MPFR_RNDN); mpfr_neg (x, x, MPFR_RNDN); mpfr_rint (y, x, MPFR_RNDN); MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_NEG(y)); /* coverage test */ mpfr_set_prec (x, 2); mpfr_set_ui (x, 1, MPFR_RNDN); mpfr_mul_2exp (x, x, mp_bits_per_limb, MPFR_RNDN); mpfr_rint (y, x, MPFR_RNDN); MPFR_ASSERTN(mpfr_cmp (y, x) == 0); /* another coverage test */ emax = mpfr_get_emax (); set_emax (1); mpfr_set_prec (x, 3); mpfr_set_str_binary (x, "1.11E0"); mpfr_set_prec (y, 2); mpfr_rint (y, x, MPFR_RNDU); /* x rounds to 1.0E1=0.1E2 which overflows */ MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0); set_emax (emax); /* yet another */ mpfr_set_prec (x, 97); mpfr_set_prec (y, 96); mpfr_set_str_binary (x, "-0.1011111001101111000111011100011100000110110110110000000111010001000101001111101010101011010111100E97"); mpfr_rint (y, x, MPFR_RNDN); MPFR_ASSERTN(mpfr_cmp (y, x) == 0); mpfr_set_prec (x, 53); mpfr_set_prec (y, 53); mpfr_set_str_binary (x, "0.10101100000000101001010101111111000000011111010000010E-1"); mpfr_rint (y, x, MPFR_RNDU); MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0); mpfr_rint (y, x, MPFR_RNDD); MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y)); mpfr_set_prec (x, 36); mpfr_set_prec (y, 2); mpfr_set_str_binary (x, "-11000110101010111111110111001.0000100"); mpfr_rint (y, x, MPFR_RNDN); mpfr_set_str_binary (x, "-11E27"); MPFR_ASSERTN(mpfr_cmp (y, x) == 0); mpfr_set_prec (x, 39); mpfr_set_prec (y, 29); mpfr_set_str_binary (x, "-0.100010110100011010001111001001001100111E39"); mpfr_rint (y, x, MPFR_RNDN); mpfr_set_str_binary (x, "-0.10001011010001101000111100101E39"); MPFR_ASSERTN(mpfr_cmp (y, x) == 0); mpfr_set_prec (x, 46); mpfr_set_prec (y, 32); mpfr_set_str_binary (x, "-0.1011100110100101000001011111101011001001101001E32"); mpfr_rint (y, x, MPFR_RNDN); mpfr_set_str_binary (x, "-0.10111001101001010000010111111011E32"); MPFR_ASSERTN(mpfr_cmp (y, x) == 0); /* coverage test for mpfr_round */ mpfr_set_prec (x, 3); mpfr_set_str_binary (x, "1.01E1"); /* 2.5 */ mpfr_set_prec (y, 2); mpfr_round (y, x); /* since mpfr_round breaks ties away, should give 3 and not 2 as with the "round to even" rule */ MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0); /* same test for the function */ (mpfr_round) (y, x); MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0); mpfr_set_prec (x, 6); mpfr_set_prec (y, 3); mpfr_set_str_binary (x, "110.111"); mpfr_round (y, x); if (mpfr_cmp_ui (y, 7)) { printf ("Error in round(110.111)\n"); exit (1); } /* Bug found by Mark J Watkins */ mpfr_set_prec (x, 84); mpfr_set_str_binary (x, "0.110011010010001000000111101101001111111100101110010000000000000" \ "000000000000000000000E32"); mpfr_round (x, x); if (mpfr_cmp_str (x, "0.1100110100100010000001111011010100000000000000" \ "00000000000000000000000000000000000000E32", 2, MPFR_RNDN)) { printf ("Rounding error when dest=src\n"); exit (1); } mpfr_clear (x); mpfr_clear (y); } #define BASIC_TEST(F,J) \ do \ { \ int red; \ for (red = 0; red <= 1; red++) \ { \ int inex1, inex2; \ unsigned int ex_flags, flags; \ \ if (red) \ { \ set_emin (e); \ set_emax (e); \ } \ \ mpfr_clear_flags (); \ inex1 = mpfr_set_si (y, J, (mpfr_rnd_t) r); \ ex_flags = __gmpfr_flags; \ mpfr_clear_flags (); \ inex2 = mpfr_rint_##F (z, x, (mpfr_rnd_t) r); \ flags = __gmpfr_flags; \ if (! (mpfr_equal_p (y, z) && \ SAME_SIGN (inex1, inex2) && \ flags == ex_flags)) \ { \ printf ("Basic test failed on mpfr_rint_" #F \ ", prec = %d, i = %d, %s\n", prec, s * i, \ mpfr_print_rnd_mode ((mpfr_rnd_t) r)); \ printf ("i.e. x = "); \ mpfr_dump (x); \ if (red) \ printf ("with emin = emax = %d\n", e); \ printf ("Expected "); \ mpfr_dump (y); \ printf ("with inex = %d (or equivalent)\n", inex1); \ printf (" flags:"); \ flags_out (ex_flags); \ printf ("Got "); \ mpfr_dump (z); \ printf ("with inex = %d (or equivalent)\n", inex2); \ printf (" flags:"); \ flags_out (flags); \ exit (1); \ } \ } \ set_emin (emin); \ set_emax (emax); \ } \ while (0) #define BASIC_TEST2(F,J,INEX) \ do \ { \ int red; \ for (red = 0; red <= 1; red++) \ { \ int inex; \ unsigned int ex_flags, flags; \ \ if (red) \ { \ set_emin (e); \ set_emax (e); \ } \ \ mpfr_clear_flags (); \ inex = mpfr_set_si (y, J, MPFR_RNDN); \ MPFR_ASSERTN (inex == 0 || mpfr_overflow_p ()); \ ex_flags = __gmpfr_flags; \ mpfr_clear_flags (); \ inex = mpfr_##F (z, x); \ if (inex != 0) \ ex_flags |= MPFR_FLAGS_INEXACT; \ flags = __gmpfr_flags; \ if (! (mpfr_equal_p (y, z) && \ inex == (INEX) && \ flags == ex_flags)) \ { \ printf ("Basic test failed on mpfr_" #F \ ", prec = %d, i = %d\n", prec, s * i); \ printf ("i.e. x = "); \ mpfr_dump (x); \ if (red) \ printf ("with emin = emax = %d\n", e); \ printf ("Expected "); \ mpfr_dump (y); \ printf ("with inex = %d\n", (INEX)); \ printf (" flags:"); \ flags_out (ex_flags); \ printf ("Got "); \ mpfr_dump (z); \ printf ("with inex = %d\n", inex); \ printf (" flags:"); \ flags_out (flags); \ exit (1); \ } \ } \ set_emin (emin); \ set_emax (emax); \ } \ while (0) /* Test mpfr_rint_* on i/4 with |i| between 1 and 72. */ static void basic_tests (void) { mpfr_t x, y, z; int prec, s, i, r; mpfr_exp_t emin, emax; emin = mpfr_get_emin (); emax = mpfr_get_emax (); mpfr_init2 (x, 16); for (prec = MPFR_PREC_MIN; prec <= 7; prec++) { mpfr_inits2 (prec, y, z, (mpfr_ptr) 0); for (s = 1; s >= -1; s -= 2) for (i = 1; i <= 72; i++) { 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); e = mpfr_get_exp (x); RND_LOOP_NO_RNDF (r) { BASIC_TEST (trunc, s * (i/4)); 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); } mpfr_clear (x); } #if __MPFR_STDC (199901L) static void test_fct (double (*f)(double), int (*g)(), char *s, mpfr_rnd_t r) { double d, y; mpfr_t dd, yy; mpfr_init2 (dd, 53); mpfr_init2 (yy, 53); for (d = -5.0; d <= 5.0; d += 0.25) { mpfr_set_d (dd, d, r); y = (*f)(d); if (g == &mpfr_rint) mpfr_rint (yy, dd, r); else (*g)(yy, dd); mpfr_set_d (dd, y, r); if (mpfr_cmp (yy, dd)) { printf ("test_against_libc: incorrect result for %s, rnd = %s," " d = %g\ngot ", s, mpfr_print_rnd_mode (r), d); mpfr_out_str (stdout, 10, 0, yy, MPFR_RNDN); printf (" instead of %g\n", y); exit (1); } } mpfr_clear (dd); mpfr_clear (yy); } #define TEST_FCT(F) test_fct (&F, &mpfr_##F, #F, r) static void test_against_libc (void) { mpfr_rnd_t r = MPFR_RNDN; (void) r; /* avoid a warning by using r */ #if HAVE_ROUND TEST_FCT (round); #endif #if HAVE_TRUNC TEST_FCT (trunc); #endif #if HAVE_FLOOR TEST_FCT (floor); #endif #if HAVE_CEIL TEST_FCT (ceil); #endif #if HAVE_NEARBYINT for (r = 0; r < MPFR_RND_MAX ; r++) if (mpfr_set_machine_rnd_mode (r) == 0) test_fct (&nearbyint, &mpfr_rint, "rint", r); #endif } #endif static void err (const char *str, mp_size_t s, mpfr_t x, mpfr_t y, mpfr_prec_t p, mpfr_rnd_t r, int trint, int inexact) { printf ("Error: %s\ns = %u, p = %u, r = %s, trint = %d, inexact = %d\nx = ", str, (unsigned int) s, (unsigned int) p, mpfr_print_rnd_mode (r), trint, inexact); mpfr_dump (x); printf ("y = "); mpfr_dump (y); exit (1); } static void coverage_03032011 (void) { mpfr_t in, out, cmp; int status; int precIn; char strData[(GMP_NUMB_BITS * 4)+256]; precIn = GMP_NUMB_BITS * 4; mpfr_init2 (in, precIn); mpfr_init2 (out, GMP_NUMB_BITS); mpfr_init2 (cmp, GMP_NUMB_BITS); /* cmp = "0.1EprecIn+2" */ /* The buffer size is sufficient, as precIn is small in practice. */ sprintf (strData, "0.1E%d", precIn+2); mpfr_set_str_binary (cmp, strData); /* in = "0.10...01EprecIn+2" use all (precIn) significand bits */ memset ((void *)strData, '0', precIn+2); strData[1] = '.'; strData[2] = '1'; sprintf (&strData[precIn+1], "1E%d", precIn+2); mpfr_set_str_binary (in, strData); status = mpfr_rint (out, in, MPFR_RNDN); if ((mpfr_cmp (out, cmp) != 0) || (status >= 0)) { printf("mpfr_rint error :\n status is %d instead of 0\n", status); printf(" out value is "); mpfr_dump(out); printf(" instead of "); mpfr_dump(cmp); exit (1); } mpfr_clear (cmp); mpfr_clear (out); mpfr_init2 (out, GMP_NUMB_BITS); mpfr_init2 (cmp, GMP_NUMB_BITS); /* cmp = "0.10...01EprecIn+2" use all (GMP_NUMB_BITS) significand bits */ strcpy (&strData[GMP_NUMB_BITS+1], &strData[precIn+1]); mpfr_set_str_binary (cmp, strData); (MPFR_MANT(in))[2] = MPFR_LIMB_HIGHBIT; status = mpfr_rint (out, in, MPFR_RNDN); if ((mpfr_cmp (out, cmp) != 0) || (status <= 0)) { printf("mpfr_rint error :\n status is %d instead of 0\n", status); printf(" out value is\n"); mpfr_dump(out); printf(" instead of\n"); mpfr_dump(cmp); exit (1); } mpfr_clear (cmp); mpfr_clear (out); mpfr_clear (in); } #define TEST_FUNCTION mpfr_rint_trunc #define TEST_RANDOM_EMIN -20 #define TEST_RANDOM_ALWAYS_SCALE 1 #define test_generic test_generic_trunc #include "tgeneric.c" #define TEST_FUNCTION mpfr_rint_floor #define TEST_RANDOM_EMIN -20 #define TEST_RANDOM_ALWAYS_SCALE 1 #define test_generic test_generic_floor #include "tgeneric.c" #define TEST_FUNCTION mpfr_rint_ceil #define TEST_RANDOM_EMIN -20 #define TEST_RANDOM_ALWAYS_SCALE 1 #define test_generic test_generic_ceil #include "tgeneric.c" #define TEST_FUNCTION mpfr_rint_round #define TEST_RANDOM_EMIN -20 #define TEST_RANDOM_ALWAYS_SCALE 1 #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[]) { mp_size_t s; mpz_t z; mpfr_prec_t p; mpfr_t x, y, t, u, v; int r; int inexact, sign_t; tests_start_mpfr (); mpfr_init (x); mpfr_init (y); mpz_init (z); mpfr_init (t); mpfr_init (u); mpfr_init (v); mpz_set_ui (z, 1); /* the code below works for 1 <= MPFR_PREC_MIN <= 2 */ MPFR_ASSERTN(1 <= MPFR_PREC_MIN && MPFR_PREC_MIN <= 2); for (s = MPFR_PREC_MIN; s < 100; s++) { if (s > 1) { mpz_mul_2exp (z, z, 1); if (randlimb () % 2) mpz_add_ui (z, z, 1); } /* now 2^(s-1) <= z < 2^s */ mpfr_set_prec (x, s); mpfr_set_prec (t, s); mpfr_set_prec (u, s); if (mpfr_set_z (x, z, MPFR_RNDN)) { gmp_printf ("Error: mpfr_set_z should be exact (z = %Zd, s = %u)\n", z, (unsigned int) s); exit (1); } if (randlimb () % 2) mpfr_neg (x, x, MPFR_RNDN); if (randlimb () % 2) mpfr_div_2ui (x, x, randlimb () % s, MPFR_RNDN); for (p = MPFR_PREC_MIN; p < 100; p++) { int trint; mpfr_set_prec (y, p); mpfr_set_prec (v, p); for (r = 0; r < MPFR_RND_MAX ; r++) 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 = (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 if (r == MPFR_RNDD) inexact = (trint ? mpfr_floor (y, x) : mpfr_rint_floor (y, x, MPFR_RNDD)); else { MPFR_ASSERTN (r == MPFR_RNDA || r == MPFR_RNDF); 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); 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); if (inexact == 0) continue; /* end of the test for exact results */ if (((r == MPFR_RNDD || (r == MPFR_RNDZ && MPFR_IS_POS (x))) && 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); 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); } 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); } 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); cmp = mpfr_cmp_abs (t, u); if (cmp > 0) err ("faithful rounding, but not the nearest integer", s, x, y, p, (mpfr_rnd_t) r, trint, inexact); if (cmp < 0) continue; /* |t| = |u|: x is the middle of two consecutive representable integers. */ if (trint == 2) { /* halfway case for mpfr_rint in MPFR_RNDN rounding 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); if (p > 1) { /* For p > 1, if floor(x) and ceil(x) aren't both representable integers, the significand must be even. */ mpfr_sub (v, v, y, MPFR_RNDN); mpfr_abs (v, v, MPFR_RNDN); if (mpfr_cmp_ui (v, 1) != 0) { 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); } } } 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); } } } } } mpfr_clear (x); mpfr_clear (y); mpz_clear (z); mpfr_clear (t); mpfr_clear (u); mpfr_clear (v); special (); basic_tests (); coverage_03032011 (); test_generic_trunc (MPFR_PREC_MIN, 300, 20); test_generic_floor (MPFR_PREC_MIN, 300, 20); test_generic_ceil (MPFR_PREC_MIN, 300, 20); test_generic_round (MPFR_PREC_MIN, 300, 20); test_generic_roundeven (MPFR_PREC_MIN, 300, 20); #if __MPFR_STDC (199901L) if (argc > 1 && strcmp (argv[1], "-s") == 0) test_against_libc (); #endif tests_end_mpfr (); return 0; }