1 /* Test file for mpfr_{mul,div}_2{ui,si}. 2 3 Copyright 1999, 2001-2004, 2006-2020 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba 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 https://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 "mpfr-test.h" 24 25 static const char * const val[] = { 26 "1.0001@100","4.0004000000000@102", "4.0004000000000@97", 27 "1.ABF012345@-100","6.afc048d140000@-98","6.afc048d140000@-103", 28 "F.FFFFFFFFF@10000","3.fffffffffc000@10003","3.fffffffffc000@9998", 29 "1.23456789ABCDEF@42","4.8d159e26af37c@44","4.8d159e26af37c@39", 30 "17@42","5.c000000000000@45","5.c000000000000@40", 31 "42@-17","1.0800000000000@-13","1.0800000000000@-18" 32 }; 33 34 static int 35 test_mul (int i, int div, mpfr_ptr y, mpfr_srcptr x, 36 unsigned long int n, mpfr_rnd_t r) 37 { 38 return 39 i == 0 ? (div ? mpfr_div_2ui : mpfr_mul_2ui) (y, x, n, r) : 40 i == 1 ? (div ? mpfr_div_2si : mpfr_mul_2si) (y, x, n, r) : 41 i == 2 ? (div ? mpfr_mul_2si : mpfr_div_2si) (y, x, -n, r) : 42 (exit (1), 0); 43 } 44 45 static void 46 underflow (mpfr_exp_t e) 47 { 48 mpfr_t x, y, z1, z2; 49 mpfr_exp_t emin; 50 int i, k, s; 51 int prec; 52 int rnd; 53 int div; 54 int inex1, inex2; 55 unsigned int flags1, flags2; 56 57 /* Test mul_2si(x, e - k), div_2si(x, k - e) and div_2ui(x, k - e) with 58 * emin = e, x = s * (1 + i/16), i in { -1, 0, 1 }, s in { -1, 1 }, and 59 * k = 1 to 4, by comparing the result with the one of a simple division. 60 */ 61 emin = mpfr_get_emin (); 62 set_emin (e); 63 mpfr_inits2 (8, x, y, (mpfr_ptr) 0); 64 for (i = 15; i <= 17; i++) 65 for (s = 1; s >= -1; s -= 2) 66 { 67 inex1 = mpfr_set_si_2exp (x, s * i, -4, MPFR_RNDN); 68 MPFR_ASSERTN (inex1 == 0); 69 for (prec = 6; prec >= 3; prec -= 3) 70 { 71 mpfr_inits2 (prec, z1, z2, (mpfr_ptr) 0); 72 RND_LOOP_NO_RNDF (rnd) 73 for (k = 1; k <= 4; k++) 74 { 75 /* The following one is assumed to be correct. */ 76 inex1 = mpfr_mul_2si (y, x, e, MPFR_RNDN); 77 MPFR_ASSERTN (inex1 == 0); 78 inex1 = mpfr_set_ui (z1, 1 << k, MPFR_RNDN); 79 MPFR_ASSERTN (inex1 == 0); 80 mpfr_clear_flags (); 81 /* Do not use mpfr_div_ui to avoid the optimization 82 by mpfr_div_2si. */ 83 inex1 = mpfr_div (z1, y, z1, (mpfr_rnd_t) rnd); 84 flags1 = __gmpfr_flags; 85 86 for (div = 0; div <= 2; div++) 87 { 88 mpfr_clear_flags (); 89 inex2 = 90 div == 0 ? 91 mpfr_mul_2si (z2, x, e - k, (mpfr_rnd_t) rnd) : 92 div == 1 ? 93 mpfr_div_2si (z2, x, k - e, (mpfr_rnd_t) rnd) : 94 mpfr_div_2ui (z2, x, k - e, (mpfr_rnd_t) rnd); 95 flags2 = __gmpfr_flags; 96 if (flags1 == flags2 && SAME_SIGN (inex1, inex2) && 97 mpfr_equal_p (z1, z2)) 98 continue; 99 printf ("Error in underflow("); 100 if (e == MPFR_EMIN_MIN) 101 printf ("MPFR_EMIN_MIN"); 102 else if (e == emin) 103 printf ("default emin"); 104 else 105 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e); 106 printf (") with mpfr_%s,\nx = %d/16, prec = %d, k = %d," 107 " %s\n", div == 0 ? "mul_2si" : div == 1 ? 108 "div_2si" : "div_2ui", s * i, prec, k, 109 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 110 printf ("Expected "); 111 mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN); 112 printf (", inex = %d, flags = %u\n", 113 VSIGN (inex1), flags1); 114 printf ("Got "); 115 mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN); 116 printf (", inex = %d, flags = %u\n", 117 VSIGN (inex2), flags2); 118 exit (1); 119 } /* div */ 120 } /* k */ 121 mpfr_clears (z1, z2, (mpfr_ptr) 0); 122 } /* prec */ 123 } /* i */ 124 mpfr_clears (x, y, (mpfr_ptr) 0); 125 set_emin (emin); 126 } 127 128 static void 129 underflow0 (void) 130 { 131 underflow (-256); 132 if (mpfr_get_emin () != MPFR_EMIN_MIN) 133 underflow (mpfr_get_emin ()); 134 underflow (MPFR_EMIN_MIN); 135 } 136 137 static void 138 large (mpfr_exp_t e) 139 { 140 mpfr_t x, y, z; 141 mpfr_exp_t emax; 142 int inex; 143 unsigned int flags; 144 145 emax = mpfr_get_emax (); 146 set_emax (e); 147 mpfr_init2 (x, 8); 148 mpfr_init2 (y, 8); 149 mpfr_init2 (z, 4); 150 151 mpfr_set_inf (x, 1); 152 mpfr_nextbelow (x); 153 154 mpfr_mul_2si (y, x, -1, MPFR_RNDU); 155 mpfr_prec_round (y, 4, MPFR_RNDU); 156 157 mpfr_clear_flags (); 158 inex = mpfr_mul_2si (z, x, -1, MPFR_RNDU); 159 flags = __gmpfr_flags; 160 161 if (inex <= 0 || flags != MPFR_FLAGS_INEXACT || ! mpfr_equal_p (y, z)) 162 { 163 printf ("Error in large("); 164 if (e == MPFR_EMAX_MAX) 165 printf ("MPFR_EMAX_MAX"); 166 else if (e == emax) 167 printf ("default emax"); 168 else 169 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e); 170 printf (") for mpfr_mul_2si\n"); 171 printf ("Expected inex > 0, flags = %u,\n y = ", 172 (unsigned int) MPFR_FLAGS_INEXACT); 173 mpfr_dump (y); 174 printf ("Got inex = %d, flags = %u,\n y = ", 175 inex, flags); 176 mpfr_dump (z); 177 exit (1); 178 } 179 180 mpfr_clear_flags (); 181 inex = mpfr_div_2si (z, x, 1, MPFR_RNDU); 182 flags = __gmpfr_flags; 183 184 if (inex <= 0 || flags != MPFR_FLAGS_INEXACT || ! mpfr_equal_p (y, z)) 185 { 186 printf ("Error in large("); 187 if (e == MPFR_EMAX_MAX) 188 printf ("MPFR_EMAX_MAX"); 189 else if (e == emax) 190 printf ("default emax"); 191 else 192 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e); 193 printf (") for mpfr_div_2si\n"); 194 printf ("Expected inex > 0, flags = %u,\n y = ", 195 (unsigned int) MPFR_FLAGS_INEXACT); 196 mpfr_dump (y); 197 printf ("Got inex = %d, flags = %u,\n y = ", 198 inex, flags); 199 mpfr_dump (z); 200 exit (1); 201 } 202 203 mpfr_clear_flags (); 204 inex = mpfr_div_2ui (z, x, 1, MPFR_RNDU); 205 flags = __gmpfr_flags; 206 207 if (inex <= 0 || flags != MPFR_FLAGS_INEXACT || ! mpfr_equal_p (y, z)) 208 { 209 printf ("Error in large("); 210 if (e == MPFR_EMAX_MAX) 211 printf ("MPFR_EMAX_MAX"); 212 else if (e == emax) 213 printf ("default emax"); 214 else 215 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e); 216 printf (") for mpfr_div_2ui\n"); 217 printf ("Expected inex > 0, flags = %u,\n y = ", 218 (unsigned int) MPFR_FLAGS_INEXACT); 219 mpfr_dump (y); 220 printf ("Got inex = %d, flags = %u,\n y = ", 221 inex, flags); 222 mpfr_dump (z); 223 exit (1); 224 } 225 226 mpfr_clears (x, y, z, (mpfr_ptr) 0); 227 set_emax (emax); 228 } 229 230 static void 231 large0 (void) 232 { 233 mpfr_exp_t emin; 234 235 emin = mpfr_get_emin (); 236 237 while (1) 238 { 239 large (256); 240 if (mpfr_get_emax () != MPFR_EMAX_MAX) 241 large (mpfr_get_emax ()); 242 large (MPFR_EMAX_MAX); 243 if (mpfr_get_emin () == MPFR_EMIN_MIN) 244 break; 245 /* Redo the test with __gmpfr_emin set to MPFR_EMIN_MIN, which can 246 be useful to trigger integer overflows as in div_2ui.c r12272. */ 247 set_emin (MPFR_EMIN_MIN); 248 } 249 250 set_emin (emin); 251 } 252 253 /* Cases where the function overflows on n = 0 when rounding is like 254 away from zero. */ 255 static void 256 overflow0 (mpfr_exp_t emax) 257 { 258 mpfr_exp_t old_emax; 259 mpfr_t x, y1, y2; 260 int neg, r, op; 261 static const char *sop[4] = { "mul_2ui", "mul_2si", "div_2ui", "div_2si" }; 262 263 old_emax = mpfr_get_emax (); 264 set_emax (emax); 265 266 mpfr_init2 (x, 8); 267 mpfr_inits2 (6, y1, y2, (mpfr_ptr) 0); 268 269 mpfr_set_inf (x, 1); 270 mpfr_nextbelow (x); 271 272 for (neg = 0; neg <= 1; neg++) 273 { 274 RND_LOOP_NO_RNDF (r) 275 { 276 int inex1, inex2; 277 mpfr_flags_t flags1, flags2; 278 279 /* Even if there isn't an overflow (rounding ~ toward zero), 280 the result is the same as the one of an overflow. */ 281 inex1 = mpfr_overflow (y1, (mpfr_rnd_t) r, neg ? -1 : 1); 282 flags1 = MPFR_FLAGS_INEXACT; 283 if (mpfr_inf_p (y1)) 284 flags1 |= MPFR_FLAGS_OVERFLOW; 285 for (op = 0; op < 4; op++) 286 { 287 mpfr_clear_flags (); 288 inex2 = 289 op == 0 ? mpfr_mul_2ui (y2, x, 0, (mpfr_rnd_t) r) : 290 op == 1 ? mpfr_mul_2si (y2, x, 0, (mpfr_rnd_t) r) : 291 op == 2 ? mpfr_div_2ui (y2, x, 0, (mpfr_rnd_t) r) : 292 op == 3 ? mpfr_div_2si (y2, x, 0, (mpfr_rnd_t) r) : 293 (MPFR_ASSERTN (0), 0); 294 flags2 = __gmpfr_flags; 295 if (!(mpfr_equal_p (y1, y2) && 296 SAME_SIGN (inex1, inex2) && 297 flags1 == flags2)) 298 { 299 printf ("Error in overflow0 for %s, mpfr_%s, emax = %" 300 MPFR_EXP_FSPEC "d,\nx = ", 301 mpfr_print_rnd_mode ((mpfr_rnd_t) r), sop[op], 302 (mpfr_eexp_t) emax); 303 mpfr_dump (x); 304 printf ("Expected "); 305 mpfr_dump (y1); 306 printf (" with inex = %d, flags =", inex1); 307 flags_out (flags1); 308 printf ("Got "); 309 mpfr_dump (y2); 310 printf (" with inex = %d, flags =", inex2); 311 flags_out (flags2); 312 exit (1); 313 } 314 } 315 } 316 mpfr_neg (x, x, MPFR_RNDN); 317 } 318 319 mpfr_clears (x, y1, y2, (mpfr_ptr) 0); 320 set_emax (old_emax); 321 } 322 323 static void 324 coverage_div_2ui (void) 325 { 326 mpfr_t x, y; 327 328 mpfr_init2 (x, 2); 329 mpfr_init2 (y, 2); 330 mpfr_set_ui_2exp (x, 1, mpfr_get_emax () - 1, MPFR_RNDN); 331 mpfr_div_2ui (y, x, (unsigned long) LONG_MAX + 1, MPFR_RNDN); 332 MPFR_ASSERTN(mpfr_zero_p (y)); 333 MPFR_ASSERTN(mpfr_signbit (y) == 0); 334 mpfr_clear (x); 335 mpfr_clear (y); 336 } 337 338 int 339 main (int argc, char *argv[]) 340 { 341 mpfr_t w,z; 342 unsigned long k; 343 int i; 344 345 tests_start_mpfr (); 346 347 coverage_div_2ui (); 348 mpfr_inits2 (53, w, z, (mpfr_ptr) 0); 349 350 for (i = 0; i < 3; i++) 351 { 352 mpfr_set_inf (w, 1); 353 test_mul (i, 0, w, w, 10, MPFR_RNDZ); 354 if (!MPFR_IS_INF(w)) 355 { 356 printf ("Result is not Inf (i = %d)\n", i); 357 exit (1); 358 } 359 360 mpfr_set_nan (w); 361 test_mul (i, 0, w, w, 10, MPFR_RNDZ); 362 if (!MPFR_IS_NAN(w)) 363 { 364 printf ("Result is not NaN (i = %d)\n", i); 365 exit (1); 366 } 367 368 for (k = 0 ; k < numberof(val) ; k+=3) 369 { 370 mpfr_set_str (w, val[k], 16, MPFR_RNDN); 371 test_mul (i, 0, z, w, 10, MPFR_RNDZ); 372 if (mpfr_cmp_str (z, val[k+1], 16, MPFR_RNDN)) 373 { 374 printf ("ERROR for x * 2^n (i = %d) for %s\n", i, val[k]); 375 printf ("Expected: %s\n" 376 "Got : ", val[k+1]); 377 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN); 378 putchar ('\n'); 379 exit (1); 380 } 381 test_mul (i, 1, z, w, 10, MPFR_RNDZ); 382 if (mpfr_cmp_str (z, val[k+2], 16, MPFR_RNDN)) 383 { 384 printf ("ERROR for x / 2^n (i = %d) for %s\n", i, val[k]); 385 printf ("Expected: %s\n" 386 "Got : ", val[k+2]); 387 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN); 388 putchar ('\n'); 389 exit (1); 390 } 391 } 392 393 mpfr_set_inf (w, 1); 394 mpfr_nextbelow (w); 395 test_mul (i, 0, w, w, 1, MPFR_RNDN); 396 if (!mpfr_inf_p (w)) 397 { 398 printf ("Overflow error (i = %d)!\n", i); 399 exit (1); 400 } 401 mpfr_set_ui (w, 0, MPFR_RNDN); 402 mpfr_nextabove (w); 403 test_mul (i, 1, w, w, 1, MPFR_RNDN); 404 if (mpfr_cmp_ui (w, 0)) 405 { 406 printf ("Underflow error (i = %d)!\n", i); 407 exit (1); 408 } 409 } 410 411 if (MPFR_EXP_MAX >= LONG_MAX/2 && MPFR_EXP_MIN <= LONG_MAX/2-LONG_MAX-1) 412 { 413 unsigned long lmp1 = (unsigned long) LONG_MAX + 1; 414 415 mpfr_set_ui (w, 1, MPFR_RNDN); 416 mpfr_mul_2ui (w, w, LONG_MAX/2, MPFR_RNDZ); 417 mpfr_div_2ui (w, w, lmp1, MPFR_RNDZ); 418 mpfr_mul_2ui (w, w, lmp1 - LONG_MAX/2, MPFR_RNDZ); 419 if (!mpfr_cmp_ui (w, 1)) 420 { 421 printf ("Underflow LONG_MAX error!\n"); 422 exit (1); 423 } 424 } 425 426 mpfr_clears (w, z, (mpfr_ptr) 0); 427 428 underflow0 (); 429 large0 (); 430 431 if (mpfr_get_emax () != MPFR_EMAX_MAX) 432 overflow0 (mpfr_get_emax ()); 433 overflow0 (MPFR_EMAX_MAX); 434 overflow0 (-1); 435 436 tests_end_mpfr (); 437 return 0; 438 } 439