1/* Test file for mpfr_{mul,div}_2{ui,si}. 2 3Copyright 1999, 2001-2004, 2006-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba 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 20https://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 "mpfr-test.h" 24 25static 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 34static int 35test_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 45static void 46underflow (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 128static void 129underflow0 (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 137static void 138large (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 230static void 231large0 (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. */ 255static void 256overflow0 (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 323static void 324coverage_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 338int 339main (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