1 /* Test file for mpfr_rint, mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round. 2 3 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire and Cacao projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 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 31 static void 32 special (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 155 static void 156 test_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 187 static void 188 test_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 213 static void 214 err (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 227 int 228 main (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