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