1/* Test file for mpfr_fmma and mpfr_fmms. 2 3Copyright 2016-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 25/* check both mpfr_fmma and mpfr_fmms */ 26static void 27random_test (mpfr_ptr a, mpfr_ptr b, mpfr_ptr c, mpfr_ptr d, mpfr_rnd_t rnd) 28{ 29 mpfr_t ref, res, ab, cd; 30 int inex_ref, inex_res; 31 mpfr_prec_t p = MPFR_PREC(a); 32 33 mpfr_init2 (res, p); 34 mpfr_init2 (ref, p); 35 mpfr_init2 (ab, mpfr_get_prec (a) + mpfr_get_prec (b)); 36 mpfr_init2 (cd, mpfr_get_prec (c) + mpfr_get_prec (d)); 37 38 /* first check fmma */ 39 inex_res = mpfr_fmma (res, a, b, c, d, rnd); 40 mpfr_mul (ab, a, b, rnd); 41 mpfr_mul (cd, c, d, rnd); 42 inex_ref = mpfr_add (ref, ab, cd, rnd); 43 if (! SAME_SIGN (inex_res, inex_ref) || 44 mpfr_nan_p (res) || mpfr_nan_p (ref) || 45 ! mpfr_equal_p (res, ref)) 46 { 47 printf ("mpfr_fmma failed for p=%ld rnd=%s\n", 48 (long int) p, mpfr_print_rnd_mode (rnd)); 49 printf ("a="); mpfr_dump (a); 50 printf ("b="); mpfr_dump (b); 51 printf ("c="); mpfr_dump (c); 52 printf ("d="); mpfr_dump (d); 53 printf ("Expected\n "); 54 mpfr_dump (ref); 55 printf (" with inex = %d\n", inex_ref); 56 printf ("Got\n "); 57 mpfr_dump (res); 58 printf (" with inex = %d\n", inex_res); 59 exit (1); 60 } 61 62 /* then check fmms */ 63 inex_res = mpfr_fmms (res, a, b, c, d, rnd); 64 mpfr_mul (ab, a, b, rnd); 65 mpfr_mul (cd, c, d, rnd); 66 inex_ref = mpfr_sub (ref, ab, cd, rnd); 67 if (! SAME_SIGN (inex_res, inex_ref) || 68 mpfr_nan_p (res) || mpfr_nan_p (ref) || 69 ! mpfr_equal_p (res, ref)) 70 { 71 printf ("mpfr_fmms failed for p=%ld rnd=%s\n", 72 (long int) p, mpfr_print_rnd_mode (rnd)); 73 printf ("a="); mpfr_dump (a); 74 printf ("b="); mpfr_dump (b); 75 printf ("c="); mpfr_dump (c); 76 printf ("d="); mpfr_dump (d); 77 printf ("Expected\n "); 78 mpfr_dump (ref); 79 printf (" with inex = %d\n", inex_ref); 80 printf ("Got\n "); 81 mpfr_dump (res); 82 printf (" with inex = %d\n", inex_res); 83 exit (1); 84 } 85 86 mpfr_clear (ab); 87 mpfr_clear (cd); 88 mpfr_clear (res); 89 mpfr_clear (ref); 90} 91 92static void 93random_tests (void) 94{ 95 mpfr_prec_t p; 96 int r; 97 mpfr_t a, b, c, d; 98 99 for (p = MPFR_PREC_MIN; p <= 4096; p++) 100 { 101 mpfr_inits2 (p, a, b, c, d, (mpfr_ptr) 0); 102 mpfr_urandomb (a, RANDS); 103 mpfr_urandomb (b, RANDS); 104 mpfr_urandomb (c, RANDS); 105 mpfr_urandomb (d, RANDS); 106 RND_LOOP_NO_RNDF (r) 107 random_test (a, b, c, d, (mpfr_rnd_t) r); 108 mpfr_clears (a, b, c, d, (mpfr_ptr) 0); 109 } 110} 111 112static void 113zero_tests (void) 114{ 115 mpfr_exp_t emin, emax; 116 mpfr_t a, b, c, d, res; 117 int i, r; 118 119 emin = mpfr_get_emin (); 120 emax = mpfr_get_emax (); 121 set_emin (MPFR_EMIN_MIN); 122 set_emax (MPFR_EMAX_MAX); 123 124 mpfr_inits2 (GMP_NUMB_BITS, a, b, c, d, (mpfr_ptr) 0); 125 for (i = 0; i <= 4; i++) 126 { 127 switch (i) 128 { 129 case 0: case 1: 130 mpfr_set_ui (a, i, MPFR_RNDN); 131 break; 132 case 2: 133 mpfr_setmax (a, MPFR_EMAX_MAX); 134 break; 135 case 3: 136 mpfr_setmin (a, MPFR_EMIN_MIN); 137 break; 138 case 4: 139 mpfr_setmin (a, MPFR_EMIN_MIN+1); 140 break; 141 } 142 RND_LOOP (r) 143 { 144 int j, inex; 145 mpfr_flags_t flags; 146 147 mpfr_set (b, a, MPFR_RNDN); 148 mpfr_set (c, a, MPFR_RNDN); 149 mpfr_neg (d, a, MPFR_RNDN); 150 /* We also want to test cases where the precision of the 151 result is twice the precision of each input, in case 152 add1sp.c and/or sub1sp.c could be involved. */ 153 for (j = 1; j <= 2; j++) 154 { 155 mpfr_init2 (res, GMP_NUMB_BITS * j); 156 mpfr_clear_flags (); 157 inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r); 158 flags = __gmpfr_flags; 159 if (flags != 0 || inex != 0 || ! mpfr_zero_p (res) || 160 (r == MPFR_RNDD ? MPFR_IS_POS (res) : MPFR_IS_NEG (res))) 161 { 162 printf ("Error in zero_tests on i = %d, %s\n", 163 i, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 164 printf ("Expected %c0, inex = 0\n", 165 r == MPFR_RNDD ? '-' : '+'); 166 printf ("Got "); 167 if (MPFR_IS_POS (res)) 168 printf ("+"); 169 mpfr_out_str (stdout, 16, 0, res, MPFR_RNDN); 170 printf (", inex = %d\n", inex); 171 printf ("Expected flags:"); 172 flags_out (0); 173 printf ("Got flags: "); 174 flags_out (flags); 175 exit (1); 176 } 177 mpfr_clear (res); 178 } /* j */ 179 } /* r */ 180 } /* i */ 181 mpfr_clears (a, b, c, d, (mpfr_ptr) 0); 182 183 set_emin (emin); 184 set_emax (emax); 185} 186 187static void 188max_tests (void) 189{ 190 mpfr_exp_t emin, emax; 191 mpfr_t x, y1, y2; 192 int r; 193 int i, inex1, inex2; 194 mpfr_flags_t flags1, flags2; 195 196 emin = mpfr_get_emin (); 197 emax = mpfr_get_emax (); 198 set_emin (MPFR_EMIN_MIN); 199 set_emax (MPFR_EMAX_MAX); 200 201 mpfr_init2 (x, GMP_NUMB_BITS); 202 mpfr_setmax (x, MPFR_EMAX_MAX); 203 flags1 = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; 204 RND_LOOP (r) 205 { 206 /* We also want to test cases where the precision of the 207 result is twice the precision of each input, in case 208 add1sp.c and/or sub1sp.c could be involved. */ 209 for (i = 1; i <= 2; i++) 210 { 211 mpfr_inits2 (GMP_NUMB_BITS * i, y1, y2, (mpfr_ptr) 0); 212 inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r); 213 mpfr_clear_flags (); 214 inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r); 215 flags2 = __gmpfr_flags; 216 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && 217 mpfr_equal_p (y1, y2))) 218 { 219 printf ("Error in max_tests for %s\n", 220 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 221 printf ("Expected "); 222 mpfr_dump (y1); 223 printf (" with inex = %d, flags =", inex1); 224 flags_out (flags1); 225 printf ("Got "); 226 mpfr_dump (y2); 227 printf (" with inex = %d, flags =", inex2); 228 flags_out (flags2); 229 exit (1); 230 } 231 mpfr_clears (y1, y2, (mpfr_ptr) 0); 232 } /* i */ 233 } /* r */ 234 mpfr_clear (x); 235 236 set_emin (emin); 237 set_emax (emax); 238} 239 240/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't. 241 * a^2 + cd where a^2 overflows and cd doesn't affect the overflow 242 * (and cd may even underflow). 243 */ 244static void 245overflow_tests (void) 246{ 247 mpfr_exp_t old_emax; 248 int i; 249 250 old_emax = mpfr_get_emax (); 251 252 for (i = 0; i < 40; i++) 253 { 254 mpfr_exp_t emax, exp_a; 255 mpfr_t a, k, c, d, z1, z2; 256 mpfr_rnd_t rnd; 257 mpfr_prec_t prec_a, prec_z; 258 int add = i & 1, inex, inex1, inex2; 259 mpfr_flags_t flags1, flags2; 260 261 /* In most cases, we do the test with the maximum exponent. */ 262 emax = i % 8 != 0 ? MPFR_EMAX_MAX : 64 + RAND_BOOL (); 263 set_emax (emax); 264 exp_a = emax/2 + 32; 265 266 rnd = RND_RAND_NO_RNDF (); 267 prec_a = 64 + randlimb () % 100; 268 prec_z = MPFR_PREC_MIN + randlimb () % 160; 269 270 mpfr_init2 (a, prec_a); 271 mpfr_urandom (a, RANDS, MPFR_RNDU); 272 mpfr_set_exp (a, exp_a); 273 274 mpfr_init2 (k, prec_a - 32); 275 mpfr_urandom (k, RANDS, MPFR_RNDU); 276 mpfr_set_exp (k, exp_a - 32); 277 278 mpfr_init2 (c, prec_a + 1); 279 inex = mpfr_add (c, a, k, MPFR_RNDN); 280 MPFR_ASSERTN (inex == 0); 281 282 mpfr_init2 (d, prec_a); 283 inex = mpfr_sub (d, a, k, MPFR_RNDN); 284 MPFR_ASSERTN (inex == 0); 285 if (add) 286 mpfr_neg (d, d, MPFR_RNDN); 287 288 mpfr_init2 (z1, prec_z); 289 mpfr_clear_flags (); 290 inex1 = mpfr_sqr (z1, k, rnd); 291 flags1 = __gmpfr_flags; 292 293 mpfr_init2 (z2, prec_z); 294 mpfr_clear_flags (); 295 inex2 = add ? 296 mpfr_fmma (z2, a, a, c, d, rnd) : 297 mpfr_fmms (z2, a, a, c, d, rnd); 298 flags2 = __gmpfr_flags; 299 300 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && 301 mpfr_equal_p (z1, z2))) 302 { 303 printf ("Error 1 in overflow_tests for %s\n", 304 mpfr_print_rnd_mode (rnd)); 305 printf ("Expected "); 306 mpfr_dump (z1); 307 printf (" with inex = %d, flags =", inex1); 308 flags_out (flags1); 309 printf ("Got "); 310 mpfr_dump (z2); 311 printf (" with inex = %d, flags =", inex2); 312 flags_out (flags2); 313 exit (1); 314 } 315 316 /* c and d such that a^2 +/- cd ~= a^2 (overflow) */ 317 mpfr_urandom (c, RANDS, MPFR_RNDU); 318 mpfr_set_exp (c, RAND_BOOL () ? exp_a - 2 : __gmpfr_emin); 319 if (RAND_BOOL ()) 320 mpfr_neg (c, c, MPFR_RNDN); 321 mpfr_urandom (d, RANDS, MPFR_RNDU); 322 mpfr_set_exp (d, RAND_BOOL () ? exp_a - 2 : __gmpfr_emin); 323 if (RAND_BOOL ()) 324 mpfr_neg (d, d, MPFR_RNDN); 325 326 mpfr_clear_flags (); 327 inex1 = mpfr_sqr (z1, a, rnd); 328 flags1 = __gmpfr_flags; 329 MPFR_ASSERTN (flags1 == (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)); 330 331 mpfr_clear_flags (); 332 inex2 = add ? 333 mpfr_fmma (z2, a, a, c, d, rnd) : 334 mpfr_fmms (z2, a, a, c, d, rnd); 335 flags2 = __gmpfr_flags; 336 337 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && 338 mpfr_equal_p (z1, z2))) 339 { 340 printf ("Error 2 in overflow_tests for %s\n", 341 mpfr_print_rnd_mode (rnd)); 342 printf ("Expected "); 343 mpfr_dump (z1); 344 printf (" with inex = %d, flags =", inex1); 345 flags_out (flags1); 346 printf ("Got "); 347 mpfr_dump (z2); 348 printf (" with inex = %d, flags =", inex2); 349 flags_out (flags2); 350 printf ("a="); mpfr_dump (a); 351 printf ("c="); mpfr_dump (c); 352 printf ("d="); mpfr_dump (d); 353 exit (1); 354 } 355 356 mpfr_clears (a, k, c, d, z1, z2, (mpfr_ptr) 0); 357 } 358 359 set_emax (old_emax); 360} 361 362/* (1/2)x + (1/2)x = x tested near underflow */ 363static void 364half_plus_half (void) 365{ 366 mpfr_exp_t emin; 367 mpfr_t h, x1, x2, y; 368 int neg, r, i, inex; 369 mpfr_flags_t flags; 370 371 emin = mpfr_get_emin (); 372 set_emin (MPFR_EMIN_MIN); 373 mpfr_inits2 (4, h, x1, (mpfr_ptr) 0); 374 mpfr_init2 (x2, GMP_NUMB_BITS); 375 mpfr_set_ui_2exp (h, 1, -1, MPFR_RNDN); 376 377 for (mpfr_setmin (x1, __gmpfr_emin); 378 MPFR_GET_EXP (x1) < __gmpfr_emin + 2; 379 mpfr_nextabove (x1)) 380 { 381 inex = mpfr_set (x2, x1, MPFR_RNDN); 382 MPFR_ASSERTN (inex == 0); 383 for (neg = 0; neg < 2; neg++) 384 { 385 RND_LOOP (r) 386 { 387 /* We also want to test cases where the precision of the 388 result is twice the precision of each input, in case 389 add1sp.c and/or sub1sp.c could be involved. */ 390 for (i = 1; i <= 2; i++) 391 { 392 mpfr_init2 (y, GMP_NUMB_BITS * i); 393 mpfr_clear_flags (); 394 inex = mpfr_fmma (y, h, x2, h, x2, (mpfr_rnd_t) r); 395 flags = __gmpfr_flags; 396 if (! (flags == 0 && inex == 0 && mpfr_equal_p (y, x2))) 397 { 398 printf ("Error in half_plus_half for %s\n", 399 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 400 printf ("Expected "); 401 mpfr_dump (x2); 402 printf (" with inex = 0, flags ="); 403 flags_out (0); 404 printf ("Got "); 405 mpfr_dump (y); 406 printf (" with inex = %d, flags =", inex); 407 flags_out (flags); 408 exit (1); 409 } 410 mpfr_clear (y); 411 } 412 } 413 mpfr_neg (x2, x2, MPFR_RNDN); 414 } 415 } 416 417 mpfr_clears (h, x1, x2, (mpfr_ptr) 0); 418 set_emin (emin); 419} 420 421/* check that result has exponent in [emin, emax] 422 (see https://sympa.inria.fr/sympa/arc/mpfr/2017-04/msg00016.html) 423 Overflow detection in sub1.c was incorrect (only for UBF cases); 424 fixed in r11414. */ 425static void 426bug20170405 (void) 427{ 428 mpfr_t x, y, z; 429 430 mpfr_inits2 (866, x, y, z, (mpfr_ptr) 0); 431 432 mpfr_set_str_binary (x, "-0.10010101110110001111000110010100011001111011110001010100010000111110001010111110100001000000011010001000010000101110000000001100001011000110010000100111001100000101110000000001001101101101010110000110100010010111011001101101010011111000101100000010001100010000011100000000011110100010111011101011000101101011110110001010011001101110011101100001111000011000000011000010101010000101001001010000111101100001000001011110011000110010001100001101101001001010000111100101000010111001001101010011001110110001000010101001100000101010110000100100100010101011111001001100010001010110011000000001011110011000110001000100101000111010010111111110010111001101110101010010101101010100111001011100101101010111010011001000001010011001010001101000111011010010100110011001111111000011101111001010111001001011011011110101101001100011010001010110011100001101100100001001100111010100010100E768635654"); 433 mpfr_set_str_binary (y, "-0.11010001010111110010110101010011000010010011010011011101100100110000110101100110111010001001110101110000011101100010110100100011001101111010100011111001011100111101110101101001000101011110101101101011010100110010111111010011011100101111110011001001010101011101111100011101100001010010011000110010110101001110010001101111111001100100000101010100110011101101101010011001000110100001001100000010110010101111000110110000111011000110001000100100100101111110001111100101011100100100110111010000010110110001110010001001101000000110111000101000110101111110000110001110100010101111010110001111010111111111010011001001100110011000110010110011000101110001010001101000100010000110011101010010010011110100000111100000101100110001111010000100011111000001101111110100000011011110010100010010011111111000010110000000011010011001100110001110111111010111110000111110010110011001000010E768635576"); 434 /* since emax = 1073741821, x^2-y^2 should overflow */ 435 mpfr_fmms (z, x, x, y, y, MPFR_RNDN); 436 MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0); 437 438 /* same test when z has a different precision */ 439 mpfr_set_prec (z, 867); 440 mpfr_fmms (z, x, x, y, y, MPFR_RNDN); 441 MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0); 442 443 mpfr_set_prec (x, 564); 444 mpfr_set_prec (y, 564); 445 mpfr_set_prec (z, 2256); 446 mpfr_set_str_binary (x, "1e-536870913"); 447 mpfr_set_str_binary (y, "-1e-536870913"); 448 mpfr_fmma (z, x, y, x, y, MPFR_RNDN); 449 /* we should get -0 as result */ 450 MPFR_ASSERTN(mpfr_zero_p (z) && mpfr_signbit (z)); 451 452 mpfr_set_prec (x, 564); 453 mpfr_set_prec (y, 564); 454 mpfr_set_prec (z, 2256); 455 mpfr_set_str_binary (x, "1.10010000111100110011001101111111101000111001011000110100110010000101000100010001000000111100010000101001011011111001111000110101111100101111001100001100011101100100011110000000011000010110111100111000100101010001011111010111011001110010001011101111001011001110110000010000011100010001010001011100100110111110101001001111001011101111110011101110101010110100010010111011111100010101111100011110111001011111101110101101101110100101111010000101011110100000000110111101000001100001000100010110100111010011011010110011100111010000101110010101111001011100110101100001100e-737194993"); 456 mpfr_set_str_binary (y, "-1.00101000100001001101011110100010110011101010011011010111100110101011111100000100101000111010111101100100110010001110011011100100110110000001011001000111101111101111110101100110111000000011000001101001010100010010001110001000011010000100111001001100101111111100010101110101001101101101111010100011011110001000010000010100011000011000010110101100000111111110111001100100100001101011111011100101110111000100101010110100010011101010110010100110100111000000100111101101101000000011110000100110100100011000010011110010001010000110100011111101101101110001110001101101010e-737194903"); 457 mpfr_fmma (z, x, y, x, y, MPFR_RNDN); 458 /* we should get -0 as result */ 459 MPFR_ASSERTN(mpfr_zero_p (z) && mpfr_signbit (z)); 460 461 mpfr_set_prec (x, 2); 462 mpfr_set_prec (y, 2); 463 mpfr_set_prec (z, 2); 464 /* (a+i*b)*(c+i*d) with: 465 a=0.10E1 466 b=0.10E-536870912 467 c=0.10E-536870912 468 d=0.10E1 */ 469 mpfr_set_str_binary (x, "0.10E1"); /* x = a = d */ 470 mpfr_set_str_binary (y, "0.10E-536870912"); /* y = b = c */ 471 /* real part is a*c-b*d = x*y-y*x */ 472 mpfr_fmms (z, x, y, y, x, MPFR_RNDN); 473 MPFR_ASSERTN(mpfr_zero_p (z) && !mpfr_signbit (z)); 474 /* imaginary part is a*d+b*c = x*x+y*y */ 475 mpfr_fmma (z, x, x, y, y, MPFR_RNDN); 476 MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0); 477 mpfr_fmma (z, y, y, x, x, MPFR_RNDN); 478 MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0); 479 480 mpfr_clears (x, y, z, (mpfr_ptr) 0); 481} 482 483static void 484underflow_tests (void) 485{ 486 mpfr_t x, y, z; 487 mpfr_prec_t p; 488 mpfr_exp_t emin; 489 490 emin = mpfr_get_emin (); 491 set_emin (-17); 492 for (p = MPFR_PREC_MIN; p <= 1024; p++) 493 { 494 mpfr_inits2 (p, x, y, (mpfr_ptr) 0); 495 mpfr_init2 (z, p + 1); 496 mpfr_set_str_binary (x, "1e-10"); 497 mpfr_set_str_binary (y, "1e-11"); 498 mpfr_clear_underflow (); 499 mpfr_fmms (z, x, x, y, y, MPFR_RNDN); 500 /* the exact result is 2^-20-2^-22, and should underflow */ 501 MPFR_ASSERTN(mpfr_underflow_p ()); 502 MPFR_ASSERTN(mpfr_zero_p (z)); 503 MPFR_ASSERTN(mpfr_signbit (z) == 0); 504 mpfr_clears (x, y, z, (mpfr_ptr) 0); 505 } 506 set_emin (emin); 507} 508 509static void 510bug20180604 (void) 511{ 512 mpfr_t x, y, yneg, z; 513 mpfr_exp_t emin; 514 int ret; 515 516 emin = mpfr_get_emin (); 517 set_emin (-1073741821); 518 mpfr_inits2 (564, x, y, yneg, (mpfr_ptr) 0); 519 mpfr_init2 (z, 2256); 520 mpfr_set_str_binary (x, "1.10010000111100110011001101111111101000111001011000110100110010000101000100010001000000111100010000101001011011111001111000110101111100101111001100001100011101100100011110000000011000010110111100111000100101010001011111010111011001110010001011101111001011001110110000010000011100010001010001011100100110111110101001001111001011101111110011101110101010110100010010111011111100010101111100011110111001011111101110101101101110100101111010000101011110100000000110111101000001100001000100010110100111010011011010110011100111010000101110010101111001011100110101100001100E-737194993"); 521 mpfr_set_str_binary (y, "-1.00101000100001001101011110100010110011101010011011010111100110101011111100000100101000111010111101100100110010001110011011100100110110000001011001000111101111101111110101100110111000000011000001101001010100010010001110001000011010000100111001001100101111111100010101110101001101101101111010100011011110001000010000010100011000011000010110101100000111111110111001100100100001101011111011100101110111000100101010110100010011101010110010100110100111000000100111101101101000000011110000100110100100011000010011110010001010000110100011111101101101110001110001101101010E-737194903"); 522 523 mpfr_clear_underflow (); 524 ret = mpfr_fmms (z, x, x, y, y, MPFR_RNDN); 525 MPFR_ASSERTN(mpfr_underflow_p ()); 526 MPFR_ASSERTN(mpfr_zero_p (z)); 527 MPFR_ASSERTN(mpfr_signbit (z) == 1); 528 MPFR_ASSERTN(ret > 0); 529 530 mpfr_clear_underflow (); 531 ret = mpfr_fmms (z, y, y, x, x, MPFR_RNDN); 532 MPFR_ASSERTN(mpfr_underflow_p ()); 533 MPFR_ASSERTN(mpfr_zero_p (z)); 534 MPFR_ASSERTN(mpfr_signbit (z) == 0); 535 MPFR_ASSERTN(ret < 0); 536 537 mpfr_neg (yneg, y, MPFR_RNDN); 538 mpfr_clear_underflow (); 539 ret = mpfr_fmms (z, x, x, y, yneg, MPFR_RNDN); 540 MPFR_ASSERTN(mpfr_underflow_p ()); 541 MPFR_ASSERTN(mpfr_zero_p (z)); 542 MPFR_ASSERTN(mpfr_signbit (z) == 0); 543 MPFR_ASSERTN(ret < 0); 544 545 mpfr_clear_underflow (); 546 ret = mpfr_fmms (z, y, yneg, x, x, MPFR_RNDN); 547 MPFR_ASSERTN(mpfr_underflow_p ()); 548 MPFR_ASSERTN(mpfr_zero_p (z)); 549 MPFR_ASSERTN(mpfr_signbit (z) == 1); 550 MPFR_ASSERTN(ret > 0); 551 552 mpfr_clears (x, y, yneg, z, (mpfr_ptr) 0); 553 set_emin (emin); 554} 555 556/* this bug was discovered from mpc_mul */ 557static void 558bug20200206 (void) 559{ 560 mpfr_exp_t emin = mpfr_get_emin (); 561 mpfr_t xre, xim, yre, yim, zre; 562 563 set_emin (-1073); 564 mpfr_inits2 (53, xre, xim, yre, yim, zre, (mpfr_ptr) 0); 565 mpfr_set_str (xre, "-0x8.294611b331c8p-904", 16, MPFR_RNDN); 566 mpfr_set_str (xim, "-0x1.859278c2992fap-676", 16, MPFR_RNDN); 567 mpfr_set_str (yre, "0x9.ac54802a95f8p-820", 16, MPFR_RNDN); 568 mpfr_set_str (yim, "0x3.17e59e7612aap-412", 16, MPFR_RNDN); 569 mpfr_fmms (zre, xre, yre, xim, yim, MPFR_RNDN); 570 if (mpfr_regular_p (zre) && mpfr_get_exp (zre) < -1073) 571 { 572 printf ("Error, mpfr_fmms returns an out-of-range exponent:\n"); 573 mpfr_dump (zre); 574 exit (1); 575 } 576 mpfr_clears (xre, xim, yre, yim, zre, (mpfr_ptr) 0); 577 set_emin (emin); 578} 579 580/* check for integer overflow (see bug fixed in r13798) */ 581static void 582extreme_underflow (void) 583{ 584 mpfr_t x, y, z; 585 mpfr_exp_t emin = mpfr_get_emin (); 586 587 set_emin (MPFR_EMIN_MIN); 588 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0); 589 mpfr_set_ui_2exp (x, 1, MPFR_EMIN_MIN - 1, MPFR_RNDN); 590 mpfr_set (y, x, MPFR_RNDN); 591 mpfr_nextabove (x); 592 mpfr_clear_flags (); 593 mpfr_fmms (z, x, x, y, y, MPFR_RNDN); 594 MPFR_ASSERTN (__gmpfr_flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT)); 595 MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z)); 596 mpfr_clears (x, y, z, (mpfr_ptr) 0); 597 set_emin (emin); 598} 599 600/* Test double-rounding cases in mpfr_set_1_2, which is called when 601 all the precisions are the same. With set.c before r13347, this 602 triggers errors for neg=0. */ 603static void 604double_rounding (void) 605{ 606 int p; 607 608 for (p = 4; p < 4 * GMP_NUMB_BITS; p++) 609 { 610 mpfr_t a, b, c, d, z1, z2; 611 int neg; 612 613 mpfr_inits2 (p, a, b, c, d, z1, z2, (mpfr_ptr) 0); 614 /* Choose a, b, c, d such that a*b = 2^p and c*d = 1 + epsilon */ 615 mpfr_set_ui (a, 2, MPFR_RNDN); 616 mpfr_set_ui_2exp (b, 1, p - 1, MPFR_RNDN); 617 mpfr_set_ui (c, 1, MPFR_RNDN); 618 mpfr_nextabove (c); 619 /* c = 1 + ulp(1) = 1 + 2 * ulp(1/2) */ 620 621 for (neg = 0; neg <= 1; neg++) 622 { 623 int inex1, inex2; 624 625 /* Set d = 1 - (1 + neg) * ulp(1/2), thus 626 * c*d = 1 + (1 - neg) * ulp(1/2) - 2 * (1 + neg) * ulp(1/2)^2, 627 * so that a*b + c*d rounds to 2^p + 1 and epsilon has the 628 * same sign as (-1)^neg. 629 */ 630 mpfr_set_ui (d, 1, MPFR_RNDN); 631 mpfr_nextbelow (d); 632 mpfr_set_ui_2exp (z1, 1, p, MPFR_RNDN); 633 if (neg) 634 { 635 mpfr_nextbelow (d); 636 inex1 = -1; 637 } 638 else 639 { 640 mpfr_nextabove (z1); 641 inex1 = 1; 642 } 643 644 inex2 = mpfr_fmma (z2, a, b, c, d, MPFR_RNDN); 645 if (!(mpfr_equal_p (z1, z2) && SAME_SIGN (inex1, inex2))) 646 { 647 printf ("Error in double_rounding for precision %d, neg=%d\n", 648 p, neg); 649 printf ("Expected "); 650 mpfr_dump (z1); 651 printf ("with inex = %d\n", inex1); 652 printf ("Got "); 653 mpfr_dump (z2); 654 printf ("with inex = %d\n", inex2); 655 exit (1); 656 } 657 } 658 659 mpfr_clears (a, b, c, d, z1, z2, (mpfr_ptr) 0); 660 } 661} 662 663int 664main (int argc, char *argv[]) 665{ 666 tests_start_mpfr (); 667 668 bug20200206 (); 669 bug20180604 (); 670 underflow_tests (); 671 random_tests (); 672 zero_tests (); 673 max_tests (); 674 overflow_tests (); 675 half_plus_half (); 676 bug20170405 (); 677 double_rounding (); 678 extreme_underflow (); 679 680 tests_end_mpfr (); 681 return 0; 682} 683