1/* Test file for mpfr_rint, mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round. 2 3Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include <stdlib.h> 24 25#include "mpfr-test.h" 26 27#if __MPFR_STDC (199901L) 28# include <math.h> 29#endif 30 31static void 32special (void) 33{ 34 mpfr_t x, y; 35 mpfr_exp_t emax; 36 37 mpfr_init (x); 38 mpfr_init (y); 39 40 mpfr_set_nan (x); 41 mpfr_rint (y, x, MPFR_RNDN); 42 MPFR_ASSERTN(mpfr_nan_p (y)); 43 44 mpfr_set_inf (x, 1); 45 mpfr_rint (y, x, MPFR_RNDN); 46 MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0); 47 48 mpfr_set_inf (x, -1); 49 mpfr_rint (y, x, MPFR_RNDN); 50 MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) < 0); 51 52 mpfr_set_ui (x, 0, MPFR_RNDN); 53 mpfr_rint (y, x, MPFR_RNDN); 54 MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y)); 55 56 mpfr_set_ui (x, 0, MPFR_RNDN); 57 mpfr_neg (x, x, MPFR_RNDN); 58 mpfr_rint (y, x, MPFR_RNDN); 59 MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_NEG(y)); 60 61 /* coverage test */ 62 mpfr_set_prec (x, 2); 63 mpfr_set_ui (x, 1, MPFR_RNDN); 64 mpfr_mul_2exp (x, x, mp_bits_per_limb, MPFR_RNDN); 65 mpfr_rint (y, x, MPFR_RNDN); 66 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 67 68 /* another coverage test */ 69 emax = mpfr_get_emax (); 70 set_emax (1); 71 mpfr_set_prec (x, 3); 72 mpfr_set_str_binary (x, "1.11E0"); 73 mpfr_set_prec (y, 2); 74 mpfr_rint (y, x, MPFR_RNDU); /* x rounds to 1.0E1=0.1E2 which overflows */ 75 MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0); 76 set_emax (emax); 77 78 /* yet another */ 79 mpfr_set_prec (x, 97); 80 mpfr_set_prec (y, 96); 81 mpfr_set_str_binary (x, "-0.1011111001101111000111011100011100000110110110110000000111010001000101001111101010101011010111100E97"); 82 mpfr_rint (y, x, MPFR_RNDN); 83 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 84 85 mpfr_set_prec (x, 53); 86 mpfr_set_prec (y, 53); 87 mpfr_set_str_binary (x, "0.10101100000000101001010101111111000000011111010000010E-1"); 88 mpfr_rint (y, x, MPFR_RNDU); 89 MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0); 90 mpfr_rint (y, x, MPFR_RNDD); 91 MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y)); 92 93 mpfr_set_prec (x, 36); 94 mpfr_set_prec (y, 2); 95 mpfr_set_str_binary (x, "-11000110101010111111110111001.0000100"); 96 mpfr_rint (y, x, MPFR_RNDN); 97 mpfr_set_str_binary (x, "-11E27"); 98 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 99 100 mpfr_set_prec (x, 39); 101 mpfr_set_prec (y, 29); 102 mpfr_set_str_binary (x, "-0.100010110100011010001111001001001100111E39"); 103 mpfr_rint (y, x, MPFR_RNDN); 104 mpfr_set_str_binary (x, "-0.10001011010001101000111100101E39"); 105 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 106 107 mpfr_set_prec (x, 46); 108 mpfr_set_prec (y, 32); 109 mpfr_set_str_binary (x, "-0.1011100110100101000001011111101011001001101001E32"); 110 mpfr_rint (y, x, MPFR_RNDN); 111 mpfr_set_str_binary (x, "-0.10111001101001010000010111111011E32"); 112 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 113 114 /* coverage test for mpfr_round */ 115 mpfr_set_prec (x, 3); 116 mpfr_set_str_binary (x, "1.01E1"); /* 2.5 */ 117 mpfr_set_prec (y, 2); 118 mpfr_round (y, x); 119 /* since mpfr_round breaks ties away, should give 3 and not 2 as with 120 the "round to even" rule */ 121 MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0); 122 /* same test for the function */ 123 (mpfr_round) (y, x); 124 MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0); 125 126 mpfr_set_prec (x, 6); 127 mpfr_set_prec (y, 3); 128 mpfr_set_str_binary (x, "110.111"); 129 mpfr_round (y, x); 130 if (mpfr_cmp_ui (y, 7)) 131 { 132 printf ("Error in round(110.111)\n"); 133 exit (1); 134 } 135 136 /* Bug found by Mark J Watkins */ 137 mpfr_set_prec (x, 84); 138 mpfr_set_str_binary (x, 139 "0.110011010010001000000111101101001111111100101110010000000000000" \ 140 "000000000000000000000E32"); 141 mpfr_round (x, x); 142 if (mpfr_cmp_str (x, "0.1100110100100010000001111011010100000000000000" \ 143 "00000000000000000000000000000000000000E32", 2, MPFR_RNDN)) 144 { 145 printf ("Rounding error when dest=src\n"); 146 exit (1); 147 } 148 149 mpfr_clear (x); 150 mpfr_clear (y); 151} 152 153#if __MPFR_STDC (199901L) 154 155static void 156test_fct (double (*f)(double), int (*g)(), char *s, mpfr_rnd_t r) 157{ 158 double d, y; 159 mpfr_t dd, yy; 160 161 mpfr_init2 (dd, 53); 162 mpfr_init2 (yy, 53); 163 for (d = -5.0; d <= 5.0; d += 0.25) 164 { 165 mpfr_set_d (dd, d, r); 166 y = (*f)(d); 167 if (g == &mpfr_rint) 168 mpfr_rint (yy, dd, r); 169 else 170 (*g)(yy, dd); 171 mpfr_set_d (dd, y, r); 172 if (mpfr_cmp (yy, dd)) 173 { 174 printf ("test_against_libc: incorrect result for %s, rnd = %s," 175 " d = %g\ngot ", s, mpfr_print_rnd_mode (r), d); 176 mpfr_out_str (stdout, 10, 0, yy, MPFR_RNDN); 177 printf (" instead of %g\n", y); 178 exit (1); 179 } 180 } 181 mpfr_clear (dd); 182 mpfr_clear (yy); 183} 184 185#define TEST_FCT(F) test_fct (&F, &mpfr_##F, #F, r) 186 187static void 188test_against_libc (void) 189{ 190 mpfr_rnd_t r = MPFR_RNDN; 191 192#if HAVE_ROUND 193 TEST_FCT (round); 194#endif 195#if HAVE_TRUNC 196 TEST_FCT (trunc); 197#endif 198#if HAVE_FLOOR 199 TEST_FCT (floor); 200#endif 201#if HAVE_CEIL 202 TEST_FCT (ceil); 203#endif 204#if HAVE_NEARBYINT 205 for (r = 0; r < MPFR_RND_MAX ; r++) 206 if (mpfr_set_machine_rnd_mode (r) == 0) 207 test_fct (&nearbyint, &mpfr_rint, "rint", r); 208#endif 209} 210 211#endif 212 213static void 214err (char *str, mp_size_t s, mpfr_t x, mpfr_t y, mpfr_prec_t p, mpfr_rnd_t r, 215 int trint, int inexact) 216{ 217 printf ("Error: %s\ns = %u, p = %u, r = %s, trint = %d, inexact = %d\nx = ", 218 str, (unsigned int) s, (unsigned int) p, mpfr_print_rnd_mode (r), 219 trint, inexact); 220 mpfr_print_binary (x); 221 printf ("\ny = "); 222 mpfr_print_binary (y); 223 printf ("\n"); 224 exit (1); 225} 226 227int 228main (int argc, char *argv[]) 229{ 230 mp_size_t s; 231 mpz_t z; 232 mpfr_prec_t p; 233 mpfr_t x, y, t, u, v; 234 int r; 235 int inexact, sign_t; 236 237 tests_start_mpfr (); 238 239 mpfr_init (x); 240 mpfr_init (y); 241 mpz_init (z); 242 mpfr_init (t); 243 mpfr_init (u); 244 mpfr_init (v); 245 mpz_set_ui (z, 1); 246 for (s = 2; s < 100; s++) 247 { 248 /* z has exactly s bits */ 249 250 mpz_mul_2exp (z, z, 1); 251 if (randlimb () % 2) 252 mpz_add_ui (z, z, 1); 253 mpfr_set_prec (x, s); 254 mpfr_set_prec (t, s); 255 mpfr_set_prec (u, s); 256 if (mpfr_set_z (x, z, MPFR_RNDN)) 257 { 258 printf ("Error: mpfr_set_z should be exact (s = %u)\n", 259 (unsigned int) s); 260 exit (1); 261 } 262 if (randlimb () % 2) 263 mpfr_neg (x, x, MPFR_RNDN); 264 if (randlimb () % 2) 265 mpfr_div_2ui (x, x, randlimb () % s, MPFR_RNDN); 266 for (p = 2; p < 100; p++) 267 { 268 int trint; 269 mpfr_set_prec (y, p); 270 mpfr_set_prec (v, p); 271 for (r = 0; r < MPFR_RND_MAX ; r++) 272 for (trint = 0; trint < 3; trint++) 273 { 274 if (trint == 2) 275 inexact = mpfr_rint (y, x, (mpfr_rnd_t) r); 276 else if (r == MPFR_RNDN) 277 inexact = mpfr_round (y, x); 278 else if (r == MPFR_RNDZ) 279 inexact = (trint ? mpfr_trunc (y, x) : 280 mpfr_rint_trunc (y, x, MPFR_RNDZ)); 281 else if (r == MPFR_RNDU) 282 inexact = (trint ? mpfr_ceil (y, x) : 283 mpfr_rint_ceil (y, x, MPFR_RNDU)); 284 else /* r = MPFR_RNDD */ 285 inexact = (trint ? mpfr_floor (y, x) : 286 mpfr_rint_floor (y, x, MPFR_RNDD)); 287 if (mpfr_sub (t, y, x, MPFR_RNDN)) 288 err ("subtraction 1 should be exact", 289 s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 290 sign_t = mpfr_cmp_ui (t, 0); 291 if (trint != 0 && 292 (((inexact == 0) && (sign_t != 0)) || 293 ((inexact < 0) && (sign_t >= 0)) || 294 ((inexact > 0) && (sign_t <= 0)))) 295 err ("wrong inexact flag", s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 296 if (inexact == 0) 297 continue; /* end of the test for exact results */ 298 299 if (((r == MPFR_RNDD || (r == MPFR_RNDZ && MPFR_SIGN (x) > 0)) 300 && inexact > 0) || 301 ((r == MPFR_RNDU || (r == MPFR_RNDZ && MPFR_SIGN (x) < 0)) 302 && inexact < 0)) 303 err ("wrong rounding direction", 304 s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 305 if (inexact < 0) 306 { 307 mpfr_add_ui (v, y, 1, MPFR_RNDU); 308 if (mpfr_cmp (v, x) <= 0) 309 err ("representable integer between x and its " 310 "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 311 } 312 else 313 { 314 mpfr_sub_ui (v, y, 1, MPFR_RNDD); 315 if (mpfr_cmp (v, x) >= 0) 316 err ("representable integer between x and its " 317 "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 318 } 319 if (r == MPFR_RNDN) 320 { 321 int cmp; 322 if (mpfr_sub (u, v, x, MPFR_RNDN)) 323 err ("subtraction 2 should be exact", 324 s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 325 cmp = mpfr_cmp_abs (t, u); 326 if (cmp > 0) 327 err ("faithful rounding, but not the nearest integer", 328 s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 329 if (cmp < 0) 330 continue; 331 /* |t| = |u|: x is the middle of two consecutive 332 representable integers. */ 333 if (trint == 2) 334 { 335 /* halfway case for mpfr_rint in MPFR_RNDN rounding 336 mode: round to an even integer or mantissa. */ 337 mpfr_div_2ui (y, y, 1, MPFR_RNDZ); 338 if (!mpfr_integer_p (y)) 339 err ("halfway case for mpfr_rint, result isn't an" 340 " even integer", s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 341 /* If floor(x) and ceil(x) aren't both representable 342 integers, the mantissa must be even. */ 343 mpfr_sub (v, v, y, MPFR_RNDN); 344 mpfr_abs (v, v, MPFR_RNDN); 345 if (mpfr_cmp_ui (v, 1) != 0) 346 { 347 mpfr_div_2si (y, y, MPFR_EXP (y) - MPFR_PREC (y) 348 + 1, MPFR_RNDN); 349 if (!mpfr_integer_p (y)) 350 err ("halfway case for mpfr_rint, mantissa isn't" 351 " even", s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 352 } 353 } 354 else 355 { /* halfway case for mpfr_round: x must have been 356 rounded away from zero. */ 357 if ((MPFR_SIGN (x) > 0 && inexact < 0) || 358 (MPFR_SIGN (x) < 0 && inexact > 0)) 359 err ("halfway case for mpfr_round, bad rounding" 360 " direction", s, x, y, p, (mpfr_rnd_t) r, trint, inexact); 361 } 362 } 363 } 364 } 365 } 366 mpfr_clear (x); 367 mpfr_clear (y); 368 mpz_clear (z); 369 mpfr_clear (t); 370 mpfr_clear (u); 371 mpfr_clear (v); 372 373 special (); 374 375#if __MPFR_STDC (199901L) 376 if (argc > 1 && strcmp (argv[1], "-s") == 0) 377 test_against_libc (); 378#endif 379 380 tests_end_mpfr (); 381 return 0; 382} 383