1 /* Test file for mpfr_fma. 2 3 Copyright 2001-2018 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 "mpfr-test.h" 24 25 /* When a * b is exact, the FMA is equivalent to the separate operations. */ 26 static void 27 test_exact (void) 28 { 29 const char *val[] = 30 { "@NaN@", "-@Inf@", "-2", "-1", "-0", "0", "1", "2", "@Inf@" }; 31 int sv = numberof (val); 32 int i, j, k; 33 int rnd; 34 mpfr_t a, b, c, r1, r2; 35 36 mpfr_inits2 (8, a, b, c, r1, r2, (mpfr_ptr) 0); 37 38 for (i = 0; i < sv; i++) 39 for (j = 0; j < sv; j++) 40 for (k = 0; k < sv; k++) 41 RND_LOOP (rnd) 42 { 43 if (mpfr_set_str (a, val[i], 10, MPFR_RNDN) || 44 mpfr_set_str (b, val[j], 10, MPFR_RNDN) || 45 mpfr_set_str (c, val[k], 10, MPFR_RNDN) || 46 mpfr_mul (r1, a, b, (mpfr_rnd_t) rnd) || 47 mpfr_add (r1, r1, c, (mpfr_rnd_t) rnd)) 48 { 49 if (rnd == MPFR_RNDF) 50 break; 51 printf ("test_exact internal error for (%d,%d,%d,%d,%s)\n", 52 i, j, k, rnd, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 53 exit (1); 54 } 55 if (mpfr_fma (r2, a, b, c, (mpfr_rnd_t) rnd)) 56 { 57 printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be exact\n", 58 i, j, k, rnd); 59 exit (1); 60 } 61 if (MPFR_IS_NAN (r1)) 62 { 63 if (MPFR_IS_NAN (r2)) 64 continue; 65 printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be NaN\n", 66 i, j, k, rnd); 67 exit (1); 68 } 69 if (! mpfr_equal_p (r1, r2) || MPFR_SIGN (r1) != MPFR_SIGN (r2)) 70 { 71 printf ("test_exact(%d,%d,%d,%d):\nexpected ", i, j, k, rnd); 72 mpfr_out_str (stdout, 10, 0, r1, MPFR_RNDN); 73 printf ("\n got "); 74 mpfr_out_str (stdout, 10, 0, r2, MPFR_RNDN); 75 printf ("\n"); 76 exit (1); 77 } 78 } 79 80 mpfr_clears (a, b, c, r1, r2, (mpfr_ptr) 0); 81 } 82 83 static void 84 test_overflow1 (void) 85 { 86 mpfr_t x, y, z, r; 87 int inex; 88 89 mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0); 90 MPFR_SET_POS (x); 91 mpfr_setmax (x, mpfr_get_emax ()); /* x = 2^emax - ulp */ 92 mpfr_set_ui (y, 2, MPFR_RNDN); /* y = 2 */ 93 mpfr_neg (z, x, MPFR_RNDN); /* z = -x = -(2^emax - ulp) */ 94 mpfr_clear_flags (); 95 /* The intermediate multiplication x * y overflows, but x * y + z = x 96 is representable. */ 97 inex = mpfr_fma (r, x, y, z, MPFR_RNDN); 98 if (inex || ! mpfr_equal_p (r, x)) 99 { 100 printf ("Error in test_overflow1\nexpected "); 101 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 102 printf (" with inex = 0\n got "); 103 mpfr_out_str (stdout, 2, 0, r, MPFR_RNDN); 104 printf (" with inex = %d\n", inex); 105 exit (1); 106 } 107 if (mpfr_overflow_p ()) 108 { 109 printf ("Error in test_overflow1: overflow flag set\n"); 110 exit (1); 111 } 112 mpfr_clears (x, y, z, r, (mpfr_ptr) 0); 113 } 114 115 static void 116 test_overflow2 (void) 117 { 118 mpfr_t x, y, z, r; 119 int i, inex, rnd, err = 0; 120 121 mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0); 122 123 MPFR_SET_POS (x); 124 mpfr_setmin (x, mpfr_get_emax ()); /* x = 0.1@emax */ 125 mpfr_set_si (y, -2, MPFR_RNDN); /* y = -2 */ 126 /* The intermediate multiplication x * y will overflow. */ 127 128 for (i = -9; i <= 9; i++) 129 RND_LOOP_NO_RNDF (rnd) 130 { 131 int inf, overflow; 132 133 inf = rnd == MPFR_RNDN || rnd == MPFR_RNDD || rnd == MPFR_RNDA; 134 overflow = inf || i <= 0; 135 136 inex = mpfr_set_si_2exp (z, i, mpfr_get_emin (), MPFR_RNDN); 137 MPFR_ASSERTN (inex == 0); 138 139 mpfr_clear_flags (); 140 /* One has: x * y = -1@emax exactly (but not representable). */ 141 inex = mpfr_fma (r, x, y, z, (mpfr_rnd_t) rnd); 142 if (overflow ^ (mpfr_overflow_p () != 0)) 143 { 144 printf ("Error in test_overflow2 (i = %d, %s): wrong overflow" 145 " flag (should be %d)\n", i, 146 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), overflow); 147 err = 1; 148 } 149 if (mpfr_nanflag_p ()) 150 { 151 printf ("Error in test_overflow2 (i = %d, %s): NaN flag should" 152 " not be set\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 153 err = 1; 154 } 155 if (mpfr_nan_p (r)) 156 { 157 printf ("Error in test_overflow2 (i = %d, %s): got NaN\n", 158 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 159 err = 1; 160 } 161 else if (MPFR_IS_POS (r)) 162 { 163 printf ("Error in test_overflow2 (i = %d, %s): wrong sign " 164 "(+ instead of -)\n", i, 165 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 166 err = 1; 167 } 168 else if (inf && ! mpfr_inf_p (r)) 169 { 170 printf ("Error in test_overflow2 (i = %d, %s): expected -Inf," 171 " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 172 mpfr_dump (r); 173 err = 1; 174 } 175 else if (!inf && (mpfr_inf_p (r) || 176 (mpfr_nextbelow (r), ! mpfr_inf_p (r)))) 177 { 178 printf ("Error in test_overflow2 (i = %d, %s): expected -MAX," 179 " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 180 mpfr_dump (r); 181 err = 1; 182 } 183 if (inf ? inex >= 0 : inex <= 0) 184 { 185 printf ("Error in test_overflow2 (i = %d, %s): wrong inexact" 186 " flag (got %d)\n", i, 187 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex); 188 err = 1; 189 } 190 191 } 192 193 if (err) 194 exit (1); 195 mpfr_clears (x, y, z, r, (mpfr_ptr) 0); 196 } 197 198 static void 199 test_underflow1 (void) 200 { 201 mpfr_t x, y, z, r; 202 int inex, signy, signz, rnd, err = 0; 203 204 mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0); 205 206 MPFR_SET_POS (x); 207 mpfr_setmin (x, mpfr_get_emin ()); /* x = 0.1@emin */ 208 209 for (signy = -1; signy <= 1; signy += 2) 210 { 211 mpfr_set_si_2exp (y, signy, -1, MPFR_RNDN); /* |y| = 1/2 */ 212 for (signz = -3; signz <= 3; signz += 2) 213 { 214 RND_LOOP (rnd) 215 { 216 mpfr_set_si (z, signz, MPFR_RNDN); 217 if (ABS (signz) != 1) 218 mpfr_setmax (z, mpfr_get_emax ()); 219 /* |z| = 1 or 2^emax - ulp */ 220 mpfr_clear_flags (); 221 inex = mpfr_fma (r, x, y, z, (mpfr_rnd_t) rnd); 222 #define STRTU1 "Error in test_underflow1 (signy = %d, signz = %d, %s)\n " 223 if (mpfr_nanflag_p ()) 224 { 225 printf (STRTU1 "NaN flag is set\n", signy, signz, 226 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 227 err = 1; 228 } 229 if (signy < 0 && MPFR_IS_LIKE_RNDD(rnd, signz)) 230 mpfr_nextbelow (z); 231 if (signy > 0 && MPFR_IS_LIKE_RNDU(rnd, signz)) 232 mpfr_nextabove (z); 233 if ((mpfr_overflow_p () != 0) ^ (mpfr_inf_p (z) != 0)) 234 { 235 printf (STRTU1 "wrong overflow flag\n", signy, signz, 236 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 237 err = 1; 238 } 239 if (mpfr_underflow_p ()) 240 { 241 printf (STRTU1 "underflow flag is set\n", signy, signz, 242 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 243 err = 1; 244 } 245 if (! mpfr_equal_p (r, z)) 246 { 247 printf (STRTU1 "got ", signy, signz, 248 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 249 mpfr_dump (r); 250 printf (" instead of "); 251 mpfr_dump (z); 252 err = 1; 253 } 254 if (inex >= 0 && (rnd == MPFR_RNDD || 255 (rnd == MPFR_RNDZ && signz > 0) || 256 (rnd == MPFR_RNDN && signy > 0))) 257 { 258 printf (STRTU1 "ternary value = %d instead of < 0\n", 259 signy, signz, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), 260 inex); 261 err = 1; 262 } 263 if (inex <= 0 && (rnd == MPFR_RNDU || 264 (rnd == MPFR_RNDZ && signz < 0) || 265 (rnd == MPFR_RNDN && signy < 0))) 266 { 267 printf (STRTU1 "ternary value = %d instead of > 0\n", 268 signy, signz, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), 269 inex); 270 err = 1; 271 } 272 } 273 } 274 } 275 276 if (err) 277 exit (1); 278 mpfr_clears (x, y, z, r, (mpfr_ptr) 0); 279 } 280 281 static void 282 test_underflow2 (void) 283 { 284 mpfr_t x, y, z, r; 285 int b, i, inex, same, err = 0; 286 287 mpfr_inits2 (32, x, y, z, r, (mpfr_ptr) 0); 288 289 mpfr_set_si_2exp (z, 1, mpfr_get_emin (), MPFR_RNDN); /* z = 2^emin */ 290 mpfr_set_si_2exp (x, 1, mpfr_get_emin (), MPFR_RNDN); /* x = 2^emin */ 291 292 for (b = 0; b <= 1; b++) 293 { 294 for (i = 15; i <= 17; i++) 295 { 296 mpfr_set_si_2exp (y, i, -4 - MPFR_PREC (z), MPFR_RNDN); 297 /* z = 1.000...00b 298 * xy = 01111 299 * or 10000 300 * or 10001 301 */ 302 mpfr_clear_flags (); 303 inex = mpfr_fma (r, x, y, z, MPFR_RNDN); 304 #define STRTU2 "Error in test_underflow2 (b = %d, i = %d)\n " 305 if (__gmpfr_flags != MPFR_FLAGS_INEXACT) 306 { 307 printf (STRTU2 "flags = %u instead of %u\n", b, i, 308 (unsigned int) __gmpfr_flags, 309 (unsigned int) MPFR_FLAGS_INEXACT); 310 err = 1; 311 } 312 same = i == 15 || (i == 16 && b == 0); 313 if (same ? (inex >= 0) : (inex <= 0)) 314 { 315 printf (STRTU2 "incorrect ternary value (%d instead of %c 0)\n", 316 b, i, inex, same ? '<' : '>'); 317 err = 1; 318 } 319 mpfr_set (y, z, MPFR_RNDN); 320 if (!same) 321 mpfr_nextabove (y); 322 if (! mpfr_equal_p (r, y)) 323 { 324 printf (STRTU2 "expected ", b, i); 325 mpfr_dump (y); 326 printf (" got "); 327 mpfr_dump (r); 328 err = 1; 329 } 330 } 331 mpfr_nextabove (z); 332 } 333 334 if (err) 335 exit (1); 336 mpfr_clears (x, y, z, r, (mpfr_ptr) 0); 337 } 338 339 static void 340 test_underflow3 (int n) 341 { 342 mpfr_t x, y, z, t1, t2; 343 int sign, k, s, rnd, inex1, inex2; 344 mpfr_exp_t e; 345 mpfr_flags_t flags1, flags2; 346 347 mpfr_inits2 (4, x, z, t1, t2, (mpfr_ptr) 0); 348 mpfr_init2 (y, 6); 349 350 e = mpfr_get_emin () - 1; 351 352 for (sign = 1; sign >= -1; sign -= 2) 353 for (k = 1; k <= 7; k++) 354 for (s = -1; s <= 1; s++) 355 { 356 mpfr_set_si_2exp (x, sign, e, MPFR_RNDN); 357 mpfr_set_si_2exp (y, 8*k+s, -6, MPFR_RNDN); 358 mpfr_neg (z, x, MPFR_RNDN); 359 /* x = sign * 2^(emin-1) 360 y = (8 * k + s) * 2^(-6) = k / 8 + s * 2^(-6) 361 z = -x = -sign * 2^(emin-1) 362 FMA(x,y,z) = sign * ((k-8) * 2^(emin-4) + s * 2^(emin-7)) exactly. 363 Note: The purpose of the s * 2^(emin-7) term is to yield 364 double rounding when scaling for k = 4, s != 0, MPFR_RNDN. */ 365 366 RND_LOOP (rnd) 367 { 368 mpfr_clear_flags (); 369 inex1 = mpfr_set_si_2exp (t1, sign * (8*k+s-64), e-6, 370 (mpfr_rnd_t) rnd); 371 flags1 = __gmpfr_flags; 372 373 mpfr_clear_flags (); 374 inex2 = mpfr_fma (t2, x, y, z, (mpfr_rnd_t) rnd); 375 flags2 = __gmpfr_flags; 376 377 if (! (mpfr_equal_p (t1, t2) && 378 SAME_SIGN (inex1, inex2) && 379 flags1 == flags2)) 380 { 381 printf ("Error in test_underflow3, n = %d, sign = %d," 382 " k = %d, s = %d, %s\n", n, sign, k, s, 383 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 384 printf ("Expected "); 385 mpfr_dump (t1); 386 printf (" with inex ~ %d, flags =", inex1); 387 flags_out (flags1); 388 printf ("Got "); 389 mpfr_dump (t2); 390 printf (" with inex ~ %d, flags =", inex2); 391 flags_out (flags2); 392 exit (1); 393 } 394 } 395 } 396 397 mpfr_clears (x, y, z, t1, t2, (mpfr_ptr) 0); 398 } 399 400 static void 401 bug20101018 (void) 402 { 403 mpfr_t x, y, z, t, u; 404 int i; 405 406 mpfr_init2 (x, 64); 407 mpfr_init2 (y, 64); 408 mpfr_init2 (z, 64); 409 mpfr_init2 (t, 64); 410 mpfr_init2 (u, 64); 411 412 mpfr_set_str (x, "0xf.fffffffffffffffp-14766", 16, MPFR_RNDN); 413 mpfr_set_str (y, "-0xf.fffffffffffffffp+317", 16, MPFR_RNDN); 414 mpfr_set_str (z, "0x8.3ffffffffffe3ffp-14443", 16, MPFR_RNDN); 415 mpfr_set_str (t, "0x8.7ffffffffffc7ffp-14444", 16, MPFR_RNDN); 416 i = mpfr_fma (u, x, y, z, MPFR_RNDN); 417 if (! mpfr_equal_p (u, t)) 418 { 419 printf ("Wrong result in bug20101018 (a)\n"); 420 printf ("Expected "); 421 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN); 422 printf ("\nGot "); 423 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN); 424 printf ("\n"); 425 exit (1); 426 } 427 if (i <= 0) 428 { 429 printf ("Wrong ternary value in bug20101018 (a)\n"); 430 printf ("Expected > 0\n"); 431 printf ("Got %d\n", i); 432 exit (1); 433 } 434 435 mpfr_set_str (x, "-0xf.fffffffffffffffp-11420", 16, MPFR_RNDN); 436 mpfr_set_str (y, "0xf.fffffffffffffffp+9863", 16, MPFR_RNDN); 437 mpfr_set_str (z, "0x8.fffff80ffffffffp-1551", 16, MPFR_RNDN); 438 mpfr_set_str (t, "0x9.fffff01ffffffffp-1552", 16, MPFR_RNDN); 439 i = mpfr_fma (u, x, y, z, MPFR_RNDN); 440 if (! mpfr_equal_p (u, t)) 441 { 442 printf ("Wrong result in bug20101018 (b)\n"); 443 printf ("Expected "); 444 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN); 445 printf ("\nGot "); 446 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN); 447 printf ("\n"); 448 exit (1); 449 } 450 if (i <= 0) 451 { 452 printf ("Wrong ternary value in bug20101018 (b)\n"); 453 printf ("Expected > 0\n"); 454 printf ("Got %d\n", i); 455 exit (1); 456 } 457 458 mpfr_set_str (x, "0xf.fffffffffffffffp-2125", 16, MPFR_RNDN); 459 mpfr_set_str (y, "-0xf.fffffffffffffffp-6000", 16, MPFR_RNDN); 460 mpfr_set_str (z, "0x8p-8119", 16, MPFR_RNDN); 461 mpfr_set_str (t, "0x8.000000000000001p-8120", 16, MPFR_RNDN); 462 i = mpfr_fma (u, x, y, z, MPFR_RNDN); 463 if (! mpfr_equal_p (u, t)) 464 { 465 printf ("Wrong result in bug20101018 (c)\n"); 466 printf ("Expected "); 467 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN); 468 printf ("\nGot "); 469 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN); 470 printf ("\n"); 471 exit (1); 472 } 473 if (i <= 0) 474 { 475 printf ("Wrong ternary value in bug20101018 (c)\n"); 476 printf ("Expected > 0\n"); 477 printf ("Got %d\n", i); 478 exit (1); 479 } 480 481 mpfr_clear (x); 482 mpfr_clear (y); 483 mpfr_clear (z); 484 mpfr_clear (t); 485 mpfr_clear (u); 486 } 487 488 /* bug found with GMP_CHECK_RANDOMIZE=1514407719 */ 489 static void 490 bug20171219 (void) 491 { 492 mpfr_t x, y, z, t, u; 493 int i; 494 495 mpfr_init2 (x, 60); 496 mpfr_init2 (y, 60); 497 mpfr_init2 (z, 60); 498 mpfr_init2 (t, 68); 499 mpfr_init2 (u, 68); 500 501 mpfr_set_ui (x, 1, MPFR_RNDN); 502 mpfr_set_ui (y, 1, MPFR_RNDN); 503 mpfr_set_ui (z, 1, MPFR_RNDN); 504 mpfr_set_ui (t, 2, MPFR_RNDN); 505 i = mpfr_fma (u, x, y, z, MPFR_RNDA); 506 if (! mpfr_equal_p (u, t)) 507 { 508 printf ("Wrong result in bug20171219 (a)\n"); 509 printf ("Expected "); 510 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN); 511 printf ("\nGot "); 512 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN); 513 printf ("\n"); 514 exit (1); 515 } 516 if (i != 0) 517 { 518 printf ("Wrong ternary value in bug20171219\n"); 519 printf ("Expected 0\n"); 520 printf ("Got %d\n", i); 521 exit (1); 522 } 523 524 mpfr_clear (x); 525 mpfr_clear (y); 526 mpfr_clear (z); 527 mpfr_clear (t); 528 mpfr_clear (u); 529 } 530 531 int 532 main (int argc, char *argv[]) 533 { 534 mpfr_t x, y, z, s; 535 mpfr_exp_t emin, emax; 536 537 tests_start_mpfr (); 538 539 emin = mpfr_get_emin (); 540 emax = mpfr_get_emax (); 541 542 bug20171219 (); 543 bug20101018 (); 544 545 mpfr_init (x); 546 mpfr_init (s); 547 mpfr_init (y); 548 mpfr_init (z); 549 550 /* check special cases */ 551 mpfr_set_prec (x, 2); 552 mpfr_set_prec (y, 2); 553 mpfr_set_prec (z, 2); 554 mpfr_set_prec (s, 2); 555 mpfr_set_str (x, "-0.75", 10, MPFR_RNDN); 556 mpfr_set_str (y, "0.5", 10, MPFR_RNDN); 557 mpfr_set_str (z, "0.375", 10, MPFR_RNDN); 558 mpfr_fma (s, x, y, z, MPFR_RNDU); /* result is 0 */ 559 if (mpfr_cmp_ui (s, 0)) 560 { 561 printf ("Error: -0.75 * 0.5 + 0.375 should be equal to 0 for prec=2\n"); 562 printf ("got instead "); 563 mpfr_dump (s); 564 exit (1); 565 } 566 567 mpfr_set_prec (x, 27); 568 mpfr_set_prec (y, 27); 569 mpfr_set_prec (z, 27); 570 mpfr_set_prec (s, 27); 571 mpfr_set_str_binary (x, "1.11111111111111111111111111e-1"); 572 mpfr_set (y, x, MPFR_RNDN); 573 mpfr_set_str_binary (z, "-1.00011110100011001011001001e-1"); 574 if (mpfr_fma (s, x, y, z, MPFR_RNDN) >= 0) 575 { 576 printf ("Wrong inexact flag for x=y=1-2^(-27)\n"); 577 exit (1); 578 } 579 580 mpfr_set_nan (x); 581 mpfr_urandomb (y, RANDS); 582 mpfr_urandomb (z, RANDS); 583 mpfr_fma (s, x, y, z, MPFR_RNDN); 584 if (!mpfr_nan_p (s)) 585 { 586 printf ("evaluation of function in x=NAN does not return NAN\n"); 587 exit (1); 588 } 589 590 mpfr_set_nan (y); 591 mpfr_urandomb (x, RANDS); 592 mpfr_urandomb (z, RANDS); 593 mpfr_fma (s, x, y, z, MPFR_RNDN); 594 if (!mpfr_nan_p(s)) 595 { 596 printf ("evaluation of function in y=NAN does not return NAN\n"); 597 exit (1); 598 } 599 600 mpfr_set_nan (z); 601 mpfr_urandomb (y, RANDS); 602 mpfr_urandomb (x, RANDS); 603 mpfr_fma (s, x, y, z, MPFR_RNDN); 604 if (!mpfr_nan_p (s)) 605 { 606 printf ("evaluation of function in z=NAN does not return NAN\n"); 607 exit (1); 608 } 609 610 mpfr_set_inf (x, 1); 611 mpfr_set_inf (y, 1); 612 mpfr_set_inf (z, 1); 613 mpfr_fma (s, x, y, z, MPFR_RNDN); 614 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0) 615 { 616 printf ("Error for (+inf) * (+inf) + (+inf)\n"); 617 exit (1); 618 } 619 620 mpfr_set_inf (x, -1); 621 mpfr_set_inf (y, -1); 622 mpfr_set_inf (z, 1); 623 mpfr_fma (s, x, y, z, MPFR_RNDN); 624 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0) 625 { 626 printf ("Error for (-inf) * (-inf) + (+inf)\n"); 627 exit (1); 628 } 629 630 mpfr_set_inf (x, 1); 631 mpfr_set_inf (y, -1); 632 mpfr_set_inf (z, -1); 633 mpfr_fma (s, x, y, z, MPFR_RNDN); 634 if (!mpfr_inf_p (s) || mpfr_sgn (s) > 0) 635 { 636 printf ("Error for (+inf) * (-inf) + (-inf)\n"); 637 exit (1); 638 } 639 640 mpfr_set_inf (x, -1); 641 mpfr_set_inf (y, 1); 642 mpfr_set_inf (z, -1); 643 mpfr_fma (s, x, y, z, MPFR_RNDN); 644 if (!mpfr_inf_p (s) || mpfr_sgn (s) > 0) 645 { 646 printf ("Error for (-inf) * (+inf) + (-inf)\n"); 647 exit (1); 648 } 649 650 mpfr_set_inf (x, 1); 651 mpfr_set_ui (y, 0, MPFR_RNDN); 652 mpfr_urandomb (z, RANDS); 653 mpfr_fma (s, x, y, z, MPFR_RNDN); 654 if (!mpfr_nan_p (s)) 655 { 656 printf ("evaluation of function in x=INF y=0 does not return NAN\n"); 657 exit (1); 658 } 659 660 mpfr_set_inf (y, 1); 661 mpfr_set_ui (x, 0, MPFR_RNDN); 662 mpfr_urandomb (z, RANDS); 663 mpfr_fma (s, x, y, z, MPFR_RNDN); 664 if (!mpfr_nan_p (s)) 665 { 666 printf ("evaluation of function in x=0 y=INF does not return NAN\n"); 667 exit (1); 668 } 669 670 mpfr_set_inf (x, 1); 671 mpfr_urandomb (y, RANDS); /* always positive */ 672 mpfr_set_inf (z, -1); 673 mpfr_fma (s, x, y, z, MPFR_RNDN); 674 if (!mpfr_nan_p (s)) 675 { 676 printf ("evaluation of function in x=INF y>0 z=-INF does not return NAN\n"); 677 exit (1); 678 } 679 680 mpfr_set_inf (y, 1); 681 mpfr_urandomb (x, RANDS); 682 mpfr_set_inf (z, -1); 683 mpfr_fma (s, x, y, z, MPFR_RNDN); 684 if (!mpfr_nan_p (s)) 685 { 686 printf ("evaluation of function in x>0 y=INF z=-INF does not return NAN\n"); 687 exit (1); 688 } 689 690 mpfr_set_inf (x, 1); 691 do mpfr_urandomb (y, RANDS); while (MPFR_IS_ZERO(y)); 692 mpfr_urandomb (z, RANDS); 693 mpfr_fma (s, x, y, z, MPFR_RNDN); 694 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0) 695 { 696 printf ("evaluation of function in x=INF does not return INF\n"); 697 exit (1); 698 } 699 700 mpfr_set_inf (y, 1); 701 do mpfr_urandomb (x, RANDS); while (MPFR_IS_ZERO(x)); 702 mpfr_urandomb (z, RANDS); 703 mpfr_fma (s, x, y, z, MPFR_RNDN); 704 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0) 705 { 706 printf ("evaluation of function at y=INF does not return INF\n"); 707 exit (1); 708 } 709 710 mpfr_set_inf (z, 1); 711 mpfr_urandomb (x, RANDS); 712 mpfr_urandomb (y, RANDS); 713 mpfr_fma (s, x, y, z, MPFR_RNDN); 714 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0) 715 { 716 printf ("evaluation of function in z=INF does not return INF\n"); 717 exit (1); 718 } 719 720 mpfr_set_ui (x, 0, MPFR_RNDN); 721 mpfr_urandomb (y, RANDS); 722 mpfr_urandomb (z, RANDS); 723 mpfr_fma (s, x, y, z, MPFR_RNDN); 724 if (! mpfr_equal_p (s, z)) 725 { 726 printf ("evaluation of function in x=0 does not return z\n"); 727 exit (1); 728 } 729 730 mpfr_set_ui (y, 0, MPFR_RNDN); 731 mpfr_urandomb (x, RANDS); 732 mpfr_urandomb (z, RANDS); 733 mpfr_fma (s, x, y, z, MPFR_RNDN); 734 if (! mpfr_equal_p (s, z)) 735 { 736 printf ("evaluation of function in y=0 does not return z\n"); 737 exit (1); 738 } 739 740 { 741 mpfr_prec_t prec; 742 mpfr_t t, slong; 743 mpfr_rnd_t rnd; 744 int inexact, compare; 745 unsigned int n; 746 747 mpfr_prec_t p0 = 2, p1 = 200; 748 unsigned int N = 200; 749 750 mpfr_init (t); 751 mpfr_init (slong); 752 753 /* generic test */ 754 for (prec = p0; prec <= p1; prec++) 755 { 756 mpfr_set_prec (x, prec); 757 mpfr_set_prec (y, prec); 758 mpfr_set_prec (z, prec); 759 mpfr_set_prec (s, prec); 760 mpfr_set_prec (t, prec); 761 762 for (n = 0; n < N; n++) 763 { 764 mpfr_urandomb (x, RANDS); 765 mpfr_urandomb (y, RANDS); 766 mpfr_urandomb (z, RANDS); 767 768 if (randlimb () % 2) 769 mpfr_neg (x, x, MPFR_RNDN); 770 if (randlimb () % 2) 771 mpfr_neg (y, y, MPFR_RNDN); 772 if (randlimb () % 2) 773 mpfr_neg (z, z, MPFR_RNDN); 774 775 rnd = RND_RAND_NO_RNDF (); 776 mpfr_set_prec (slong, 2 * prec); 777 if (mpfr_mul (slong, x, y, rnd)) 778 { 779 printf ("x*y should be exact\n"); 780 exit (1); 781 } 782 compare = mpfr_add (t, slong, z, rnd); 783 inexact = mpfr_fma (s, x, y, z, rnd); 784 if (! mpfr_equal_p (s, t)) 785 { 786 printf ("results differ for\n"); 787 printf (" x="); 788 mpfr_dump (x); 789 printf (" y="); 790 mpfr_dump (y); 791 printf (" z="); 792 mpfr_dump (z); 793 printf (" with prec=%u rnd_mode=%s\n", (unsigned int) prec, 794 mpfr_print_rnd_mode (rnd)); 795 printf ("got "); 796 mpfr_dump (s); 797 printf ("expected "); 798 mpfr_dump (t); 799 printf ("approx "); 800 mpfr_dump (slong); 801 exit (1); 802 } 803 if (! SAME_SIGN (inexact, compare)) 804 { 805 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n", 806 mpfr_print_rnd_mode (rnd), compare, inexact); 807 printf (" x="); mpfr_dump (x); 808 printf (" y="); mpfr_dump (y); 809 printf (" z="); mpfr_dump (z); 810 printf (" s="); mpfr_dump (s); 811 exit (1); 812 } 813 } 814 } 815 mpfr_clear (t); 816 mpfr_clear (slong); 817 818 } 819 mpfr_clear (x); 820 mpfr_clear (y); 821 mpfr_clear (z); 822 mpfr_clear (s); 823 824 test_exact (); 825 826 test_overflow1 (); 827 test_overflow2 (); 828 test_underflow1 (); 829 test_underflow2 (); 830 test_underflow3 (1); 831 832 set_emin (MPFR_EMIN_MIN); 833 set_emax (MPFR_EMAX_MAX); 834 test_overflow1 (); 835 test_overflow2 (); 836 test_underflow1 (); 837 test_underflow2 (); 838 test_underflow3 (2); 839 set_emin (emin); 840 set_emax (emax); 841 842 tests_end_mpfr (); 843 return 0; 844 } 845