1/* Test file for the various power functions 2 3Copyright 2008-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/* Note: some tests of the other tpow* test files could be moved there. 24 The main goal of this test file is to test _all_ the power functions 25 on special values, to make sure that they are consistent and give the 26 expected result, in particular because such special cases are handled 27 in different ways in each function. */ 28 29/* Execute with at least an argument to report all the errors found by 30 comparisons. */ 31 32#include "mpfr-test.h" 33 34/* Behavior of cmpres (called by test_others): 35 * 0: stop as soon as an error is found. 36 * 1: report all errors found by test_others. 37 * -1: the 1 is changed to this value as soon as an error has been found. 38 */ 39static int all_cmpres_errors; 40 41/* Non-zero when extended exponent range */ 42static int ext = 0; 43 44static const char *val[] = 45 { "min", "min+", "max", "@NaN@", "-@Inf@", "-4", "-3", "-2", "-1.5", 46 "-1", "-0.5", "-0", "0", "0.5", "1", "1.5", "2", "3", "4", "@Inf@" }; 47 48static void 49err (const char *s, int i, int j, int rnd, mpfr_srcptr z, int inex) 50{ 51 puts (s); 52 if (ext) 53 puts ("extended exponent range"); 54 printf ("x = %s, y = %s, %s\n", val[i], val[j], 55 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 56 printf ("z = "); 57 mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN); 58 printf ("\ninex = %d\n", inex); 59 exit (1); 60} 61 62/* Arguments: 63 * spx: non-zero if px is a stringm zero if px is a MPFR number. 64 * px: value of x (string or MPFR number). 65 * sy: value of y (string). 66 * rnd: rounding mode. 67 * z1: expected result (null pointer if unknown pure FP value). 68 * inex1: expected ternary value (if z1 is not a null pointer). 69 * z2: computed result. 70 * inex2: computed ternary value. 71 * flags1: expected flags (computed flags in __gmpfr_flags). 72 * s1, s2: strings about the context. 73 */ 74static void 75cmpres (int spx, const void *px, const char *sy, mpfr_rnd_t rnd, 76 mpfr_srcptr z1, int inex1, mpfr_srcptr z2, int inex2, 77 unsigned int flags1, const char *s1, const char *s2) 78{ 79 unsigned int flags2 = __gmpfr_flags; 80 81 if (flags1 == flags2) 82 { 83 /* Note: the test on the sign of z1 and z2 is needed 84 in case they are both zeros. */ 85 if (z1 == NULL) 86 { 87 if (MPFR_IS_PURE_FP (z2)) 88 return; 89 } 90 else if (SAME_SIGN (inex1, inex2) && SAME_VAL (z1, z2)) 91 return; 92 } 93 94 printf ("Error in %s\nwith %s%s\nx = ", s1, s2, 95 ext ? ", extended exponent range" : ""); 96 if (spx) 97 printf ("%s, ", (char *) px); 98 else 99 { 100 mpfr_out_str (stdout, 16, 0, (mpfr_ptr) px, MPFR_RNDN); 101 puts (","); 102 } 103 printf ("y = %s, %s\n", sy, mpfr_print_rnd_mode (rnd)); 104 printf ("Expected "); 105 if (z1 == NULL) 106 { 107 printf ("pure FP value, flags ="); 108 flags_out (flags1); 109 } 110 else 111 { 112 mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN); 113 printf (", inex = %d,\n flags =", VSIGN (inex1)); 114 flags_out (flags1); 115 } 116 printf ("Got "); 117 mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN); 118 printf (", inex = %d,\n flags =", VSIGN (inex2)); 119 flags_out (flags2); 120 if (all_cmpres_errors != 0) 121 all_cmpres_errors = -1; 122 else 123 exit (1); 124} 125 126static int 127is_odd (mpfr_srcptr x) 128{ 129 /* works only with the values from val[] */ 130 return mpfr_integer_p (x) && mpfr_fits_slong_p (x, MPFR_RNDN) && 131 (mpfr_get_si (x, MPFR_RNDN) & 1); 132} 133 134/* Compare the result (z1,inex1) of mpfr_pow with all flags cleared 135 with those of mpfr_pow with all flags set and of the other power 136 functions. Arguments x and y are the input values; sx and sy are 137 their string representations (sx may be null); rnd contains the 138 rounding mode; s is a string containing the function that called 139 test_others. */ 140static void 141test_others (const void *sx, const char *sy, mpfr_rnd_t rnd, 142 mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z1, 143 int inex1, unsigned int flags, const char *s) 144{ 145 mpfr_t z2; 146 int inex2; 147 int spx = sx != NULL; 148 149 if (!spx) 150 sx = x; 151 152 mpfr_init2 (z2, mpfr_get_prec (z1)); 153 154 __gmpfr_flags = MPFR_FLAGS_ALL; 155 inex2 = mpfr_pow (z2, x, y, rnd); 156 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 157 s, "mpfr_pow, flags set"); 158 159 /* If y is an integer that fits in an unsigned long and is not -0, 160 we can test mpfr_pow_ui. */ 161 if (MPFR_IS_POS (y) && mpfr_integer_p (y) && 162 mpfr_fits_ulong_p (y, MPFR_RNDN)) 163 { 164 unsigned long yy = mpfr_get_ui (y, MPFR_RNDN); 165 166 mpfr_clear_flags (); 167 inex2 = mpfr_pow_ui (z2, x, yy, rnd); 168 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 169 s, "mpfr_pow_ui, flags cleared"); 170 __gmpfr_flags = MPFR_FLAGS_ALL; 171 inex2 = mpfr_pow_ui (z2, x, yy, rnd); 172 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 173 s, "mpfr_pow_ui, flags set"); 174 175 /* If x is an integer that fits in an unsigned long and is not -0, 176 we can also test mpfr_ui_pow_ui. */ 177 if (MPFR_IS_POS (x) && mpfr_integer_p (x) && 178 mpfr_fits_ulong_p (x, MPFR_RNDN)) 179 { 180 unsigned long xx = mpfr_get_ui (x, MPFR_RNDN); 181 182 mpfr_clear_flags (); 183 inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd); 184 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 185 s, "mpfr_ui_pow_ui, flags cleared"); 186 __gmpfr_flags = MPFR_FLAGS_ALL; 187 inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd); 188 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 189 s, "mpfr_ui_pow_ui, flags set"); 190 } 191 } 192 193 /* If y is an integer but not -0 and not huge, we can test mpfr_pow_z, 194 and possibly mpfr_pow_si (and possibly mpfr_ui_div). */ 195 if (MPFR_IS_ZERO (y) ? MPFR_IS_POS (y) : 196 (mpfr_integer_p (y) && MPFR_GET_EXP (y) < 256)) 197 { 198 mpz_t yyy; 199 200 /* If y fits in a long, we can test mpfr_pow_si. */ 201 if (mpfr_fits_slong_p (y, MPFR_RNDN)) 202 { 203 long yy = mpfr_get_si (y, MPFR_RNDN); 204 205 mpfr_clear_flags (); 206 inex2 = mpfr_pow_si (z2, x, yy, rnd); 207 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 208 s, "mpfr_pow_si, flags cleared"); 209 __gmpfr_flags = MPFR_FLAGS_ALL; 210 inex2 = mpfr_pow_si (z2, x, yy, rnd); 211 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 212 s, "mpfr_pow_si, flags set"); 213 214 /* If y = -1, we can test mpfr_ui_div. */ 215 if (yy == -1) 216 { 217 mpfr_clear_flags (); 218 inex2 = mpfr_ui_div (z2, 1, x, rnd); 219 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 220 s, "mpfr_ui_div, flags cleared"); 221 __gmpfr_flags = MPFR_FLAGS_ALL; 222 inex2 = mpfr_ui_div (z2, 1, x, rnd); 223 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 224 s, "mpfr_ui_div, flags set"); 225 } 226 227 /* If y = 2, we can test mpfr_sqr. */ 228 if (yy == 2) 229 { 230 mpfr_clear_flags (); 231 inex2 = mpfr_sqr (z2, x, rnd); 232 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 233 s, "mpfr_sqr, flags cleared"); 234 __gmpfr_flags = MPFR_FLAGS_ALL; 235 inex2 = mpfr_sqr (z2, x, rnd); 236 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 237 s, "mpfr_sqr, flags set"); 238 } 239 } 240 241 /* Test mpfr_pow_z. */ 242 mpz_init (yyy); 243 mpfr_get_z (yyy, y, MPFR_RNDN); 244 mpfr_clear_flags (); 245 inex2 = mpfr_pow_z (z2, x, yyy, rnd); 246 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 247 s, "mpfr_pow_z, flags cleared"); 248 __gmpfr_flags = MPFR_FLAGS_ALL; 249 inex2 = mpfr_pow_z (z2, x, yyy, rnd); 250 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 251 s, "mpfr_pow_z, flags set"); 252 mpz_clear (yyy); 253 } 254 255 /* If y = 0.5, we can test mpfr_sqrt, except if x is -0 or -Inf (because 256 the rule for mpfr_pow on these special values is different). */ 257 if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "0.5") == 0 && 258 ! ((MPFR_IS_ZERO (x) || MPFR_IS_INF (x)) && MPFR_IS_NEG (x))) 259 { 260 mpfr_clear_flags (); 261 inex2 = mpfr_sqrt (z2, x, rnd); 262 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 263 s, "mpfr_sqrt, flags cleared"); 264 __gmpfr_flags = MPFR_FLAGS_ALL; 265 inex2 = mpfr_sqrt (z2, x, rnd); 266 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 267 s, "mpfr_sqrt, flags set"); 268 } 269 270 /* If y = -0.5, we can test mpfr_rec_sqrt, except if x = -Inf 271 (because the rule for mpfr_pow on -Inf is different). */ 272 if (MPFR_IS_PURE_FP (y) && mpfr_cmp_str1 (y, "-0.5") == 0 && 273 ! (MPFR_IS_INF (x) && MPFR_IS_NEG (x))) 274 { 275 mpfr_clear_flags (); 276 inex2 = mpfr_rec_sqrt (z2, x, rnd); 277 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 278 s, "mpfr_rec_sqrt, flags cleared"); 279 __gmpfr_flags = MPFR_FLAGS_ALL; 280 inex2 = mpfr_rec_sqrt (z2, x, rnd); 281 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 282 s, "mpfr_rec_sqrt, flags set"); 283 } 284 285 /* If x is an integer that fits in an unsigned long and is not -0, 286 we can test mpfr_ui_pow. */ 287 if (MPFR_IS_POS (x) && mpfr_integer_p (x) && 288 mpfr_fits_ulong_p (x, MPFR_RNDN)) 289 { 290 unsigned long xx = mpfr_get_ui (x, MPFR_RNDN); 291 292 mpfr_clear_flags (); 293 inex2 = mpfr_ui_pow (z2, xx, y, rnd); 294 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 295 s, "mpfr_ui_pow, flags cleared"); 296 __gmpfr_flags = MPFR_FLAGS_ALL; 297 inex2 = mpfr_ui_pow (z2, xx, y, rnd); 298 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 299 s, "mpfr_ui_pow, flags set"); 300 301 /* If x = 2, we can test mpfr_exp2. */ 302 if (xx == 2) 303 { 304 mpfr_clear_flags (); 305 inex2 = mpfr_exp2 (z2, y, rnd); 306 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 307 s, "mpfr_exp2, flags cleared"); 308 __gmpfr_flags = MPFR_FLAGS_ALL; 309 inex2 = mpfr_exp2 (z2, y, rnd); 310 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 311 s, "mpfr_exp2, flags set"); 312 } 313 314 /* If x = 10, we can test mpfr_exp10. */ 315 if (xx == 10) 316 { 317 mpfr_clear_flags (); 318 inex2 = mpfr_exp10 (z2, y, rnd); 319 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags, 320 s, "mpfr_exp10, flags cleared"); 321 __gmpfr_flags = MPFR_FLAGS_ALL; 322 inex2 = mpfr_exp10 (z2, y, rnd); 323 cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL, 324 s, "mpfr_exp10, flags set"); 325 } 326 } 327 328 mpfr_clear (z2); 329} 330 331static int 332my_setstr (mpfr_ptr t, const char *s) 333{ 334 if (strcmp (s, "min") == 0) 335 { 336 mpfr_setmin (t, mpfr_get_emin ()); 337 MPFR_SET_POS (t); 338 return 0; 339 } 340 if (strcmp (s, "min+") == 0) 341 { 342 mpfr_setmin (t, mpfr_get_emin ()); 343 MPFR_SET_POS (t); 344 mpfr_nextabove (t); 345 return 0; 346 } 347 if (strcmp (s, "max") == 0) 348 { 349 mpfr_setmax (t, mpfr_get_emax ()); 350 MPFR_SET_POS (t); 351 return 0; 352 } 353 return mpfr_set_str (t, s, 10, MPFR_RNDN); 354} 355 356static void 357tst (void) 358{ 359 int sv = numberof (val); 360 int i, j; 361 int rnd; 362 mpfr_t x, y, z, tmp; 363 364 mpfr_inits2 (53, x, y, z, tmp, (mpfr_ptr) 0); 365 366 for (i = 0; i < sv; i++) 367 for (j = 0; j < sv; j++) 368 RND_LOOP_NO_RNDF (rnd) 369 { 370 int exact, inex; 371 unsigned int flags; 372 373 if (my_setstr (x, val[i]) || my_setstr (y, val[j])) 374 { 375 printf ("internal error for (%d,%d,%d)\n", i, j, rnd); 376 exit (1); 377 } 378 mpfr_clear_flags (); 379 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd); 380 flags = __gmpfr_flags; 381 if (! MPFR_IS_NAN (z) && mpfr_nanflag_p ()) 382 err ("got NaN flag without NaN value", i, j, rnd, z, inex); 383 if (MPFR_IS_NAN (z) && ! mpfr_nanflag_p ()) 384 err ("got NaN value without NaN flag", i, j, rnd, z, inex); 385 if (inex != 0 && ! mpfr_inexflag_p ()) 386 err ("got non-zero ternary value without inexact flag", 387 i, j, rnd, z, inex); 388 if (inex == 0 && mpfr_inexflag_p ()) 389 err ("got null ternary value with inexact flag", 390 i, j, rnd, z, inex); 391 if (i >= 3 && j >= 3) 392 { 393 if (mpfr_underflow_p ()) 394 err ("got underflow", i, j, rnd, z, inex); 395 if (mpfr_overflow_p ()) 396 err ("got overflow", i, j, rnd, z, inex); 397 exact = MPFR_IS_SINGULAR (z) || 398 (mpfr_mul_2ui (tmp, z, 16, MPFR_RNDN), mpfr_integer_p (tmp)); 399 if (exact && inex != 0) 400 err ("got exact value with ternary flag different from 0", 401 i, j, rnd, z, inex); 402 if (! exact && inex == 0) 403 err ("got inexact value with ternary flag equal to 0", 404 i, j, rnd, z, inex); 405 } 406 if (MPFR_IS_ZERO (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y)) 407 { 408 if (MPFR_IS_NEG (y) && ! MPFR_IS_INF (z)) 409 err ("expected an infinity", i, j, rnd, z, inex); 410 if (MPFR_IS_POS (y) && ! MPFR_IS_ZERO (z)) 411 err ("expected a zero", i, j, rnd, z, inex); 412 if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z)) 413 err ("wrong sign", i, j, rnd, z, inex); 414 } 415 if (! MPFR_IS_NAN (x) && mpfr_cmp_si (x, -1) == 0) 416 { 417 /* x = -1 */ 418 if (! (MPFR_IS_INF (y) || mpfr_integer_p (y)) && 419 ! MPFR_IS_NAN (z)) 420 err ("expected NaN", i, j, rnd, z, inex); 421 if ((MPFR_IS_INF (y) || (mpfr_integer_p (y) && ! is_odd (y))) 422 && ! mpfr_equal_p (z, __gmpfr_one)) 423 err ("expected 1", i, j, rnd, z, inex); 424 if (is_odd (y) && 425 (MPFR_IS_NAN (z) || mpfr_cmp_si (z, -1) != 0)) 426 err ("expected -1", i, j, rnd, z, inex); 427 } 428 if ((mpfr_equal_p (x, __gmpfr_one) || MPFR_IS_ZERO (y)) && 429 ! mpfr_equal_p (z, __gmpfr_one)) 430 err ("expected 1", i, j, rnd, z, inex); 431 if (MPFR_IS_PURE_FP (x) && MPFR_IS_NEG (x) && 432 MPFR_IS_FP (y) && ! mpfr_integer_p (y) && 433 ! MPFR_IS_NAN (z)) 434 err ("expected NaN", i, j, rnd, z, inex); 435 if (MPFR_IS_INF (y) && MPFR_NOTZERO (x)) 436 { 437 int cmpabs1 = mpfr_cmpabs (x, __gmpfr_one); 438 439 if ((MPFR_IS_NEG (y) ? (cmpabs1 < 0) : (cmpabs1 > 0)) && 440 ! (MPFR_IS_POS (z) && MPFR_IS_INF (z))) 441 err ("expected +Inf", i, j, rnd, z, inex); 442 if ((MPFR_IS_NEG (y) ? (cmpabs1 > 0) : (cmpabs1 < 0)) && 443 ! (MPFR_IS_POS (z) && MPFR_IS_ZERO (z))) 444 err ("expected +0", i, j, rnd, z, inex); 445 } 446 if (MPFR_IS_INF (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y)) 447 { 448 if (MPFR_IS_POS (y) && ! MPFR_IS_INF (z)) 449 err ("expected an infinity", i, j, rnd, z, inex); 450 if (MPFR_IS_NEG (y) && ! MPFR_IS_ZERO (z)) 451 err ("expected a zero", i, j, rnd, z, inex); 452 if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z)) 453 err ("wrong sign", i, j, rnd, z, inex); 454 } 455 test_others (val[i], val[j], (mpfr_rnd_t) rnd, x, y, z, inex, flags, 456 "tst"); 457 } 458 mpfr_clears (x, y, z, tmp, (mpfr_ptr) 0); 459} 460 461static void 462underflow_up1 (void) 463{ 464 mpfr_t delta, x, y, z, z0; 465 mpfr_exp_t n; 466 int inex; 467 int rnd; 468 int i; 469 470 n = mpfr_get_emin (); 471 if (n < LONG_MIN) 472 return; 473 474 mpfr_init2 (delta, 2); 475 inex = mpfr_set_ui_2exp (delta, 1, -2, MPFR_RNDN); 476 MPFR_ASSERTN (inex == 0); 477 478 mpfr_init2 (x, 8); 479 inex = mpfr_set_ui (x, 2, MPFR_RNDN); 480 MPFR_ASSERTN (inex == 0); 481 482 mpfr_init2 (y, sizeof (long) * CHAR_BIT + 4); 483 inex = mpfr_set_si (y, n, MPFR_RNDN); 484 MPFR_ASSERTN (inex == 0); 485 486 mpfr_init2 (z0, 2); 487 mpfr_set_ui (z0, 0, MPFR_RNDN); 488 489 mpfr_init2 (z, 32); 490 491 for (i = 0; i <= 12; i++) 492 { 493 unsigned int flags = 0; 494 char sy[256]; /* larger than needed, for maintainability */ 495 496 /* Test 2^(emin - i/4). 497 * --> Underflow iff i > 4. 498 * --> Zero in MPFR_RNDN iff i >= 8. 499 */ 500 501 if (i != 0 && i != 4) 502 flags |= MPFR_FLAGS_INEXACT; 503 if (i > 4) 504 flags |= MPFR_FLAGS_UNDERFLOW; 505 506 sprintf (sy, "emin - %d/4", i); 507 508 RND_LOOP_NO_RNDF (rnd) 509 { 510 int zero; 511 512 zero = (i > 4 && (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)) || 513 (i >= 8 && rnd == MPFR_RNDN); 514 515 mpfr_clear_flags (); 516 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd); 517 cmpres (1, "2", sy, (mpfr_rnd_t) rnd, zero ? z0 : (mpfr_ptr) NULL, 518 -1, z, inex, flags, "underflow_up1", "mpfr_pow"); 519 test_others ("2", sy, (mpfr_rnd_t) rnd, x, y, z, inex, flags, 520 "underflow_up1"); 521 } 522 523 inex = mpfr_sub (y, y, delta, MPFR_RNDN); 524 MPFR_ASSERTN (inex == 0); 525 } 526 527 mpfr_clears (delta, x, y, z, z0, (mpfr_ptr) 0); 528} 529 530/* With pow.c r5497, the following test fails on a 64-bit Linux machine 531 * due to a double-rounding problem when rescaling the result: 532 * Error with underflow_up2 and extended exponent range 533 * x = 7.fffffffffffffff0@-1, 534 * y = 4611686018427387904, MPFR_RNDN 535 * Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9 536 * Got 0, inex = -1, flags = 9 537 * With pow_ui.c r5423, the following test fails on a 64-bit Linux machine 538 * as underflows and overflows are not handled correctly (the approximation 539 * error is ignored): 540 * Error with mpfr_pow_ui, flags cleared 541 * x = 7.fffffffffffffff0@-1, 542 * y = 4611686018427387904, MPFR_RNDN 543 * Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9 544 * Got 0, inex = -1, flags = 9 545 */ 546static void 547underflow_up2 (void) 548{ 549 mpfr_t x, y, z, z0, eps; 550 mpfr_exp_t n; 551 int inex; 552 int rnd; 553 554 n = 1 - mpfr_get_emin (); 555 MPFR_ASSERTN (n > 1); 556 if (n > ULONG_MAX) 557 return; 558 559 mpfr_init2 (eps, 2); 560 mpfr_set_ui_2exp (eps, 1, -1, MPFR_RNDN); /* 1/2 */ 561 mpfr_div_ui (eps, eps, n, MPFR_RNDZ); /* 1/(2n) rounded toward zero */ 562 563 mpfr_init2 (x, sizeof (unsigned long) * CHAR_BIT + 1); 564 inex = mpfr_ui_sub (x, 1, eps, MPFR_RNDN); 565 MPFR_ASSERTN (inex == 0); /* since n < 2^(size_of_long_in_bits) */ 566 inex = mpfr_div_2ui (x, x, 1, MPFR_RNDN); /* 1/2 - eps/2 exactly */ 567 MPFR_ASSERTN (inex == 0); 568 569 mpfr_init2 (y, sizeof (unsigned long) * CHAR_BIT); 570 inex = mpfr_set_ui (y, n, MPFR_RNDN); 571 MPFR_ASSERTN (inex == 0); 572 573 /* 0 < eps < 1 / (2n), thus (1 - eps)^n > 1/2, 574 and 1/2 (1/2)^n < (1/2 - eps/2)^n < (1/2)^n. */ 575 mpfr_inits2 (64, z, z0, (mpfr_ptr) 0); 576 RND_LOOP_NO_RNDF (rnd) 577 { 578 unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; 579 int expected_inex; 580 char sy[256]; 581 582 mpfr_set_ui (z0, 0, MPFR_RNDN); 583 expected_inex = rnd == MPFR_RNDN || rnd == MPFR_RNDU || rnd == MPFR_RNDA ? 584 (mpfr_nextabove (z0), 1) : -1; 585 sprintf (sy, "%lu", (unsigned long) n); 586 587 mpfr_clear_flags (); 588 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd); 589 cmpres (0, x, sy, (mpfr_rnd_t) rnd, z0, expected_inex, z, inex, ufinex, 590 "underflow_up2", "mpfr_pow"); 591 test_others (NULL, sy, (mpfr_rnd_t) rnd, x, y, z, inex, ufinex, 592 "underflow_up2"); 593 } 594 595 mpfr_clears (x, y, z, z0, eps, (mpfr_ptr) 0); 596} 597 598static void 599underflow_up3 (void) 600{ 601 mpfr_t x, y, z, z0; 602 int inex; 603 int rnd; 604 int i; 605 606 mpfr_init2 (x, 64); 607 mpfr_init2 (y, sizeof (mpfr_exp_t) * CHAR_BIT); 608 mpfr_init2 (z, 32); 609 mpfr_init2 (z0, 2); 610 611 inex = mpfr_set_exp_t (y, mpfr_get_emin () - 2, MPFR_RNDN); 612 MPFR_ASSERTN (inex == 0); 613 for (i = -1; i <= 1; i++) 614 RND_LOOP_NO_RNDF (rnd) 615 { 616 unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; 617 int expected_inex; 618 619 mpfr_set_ui (x, 2, MPFR_RNDN); 620 if (i < 0) 621 mpfr_nextbelow (x); 622 if (i > 0) 623 mpfr_nextabove (x); 624 /* x = 2 + i * eps, y = emin - 2, x^y ~= 2^(emin - 2) */ 625 626 expected_inex = rnd == MPFR_RNDU || rnd == MPFR_RNDA 627 || (rnd == MPFR_RNDN && i < 0) ? 1 : -1; 628 629 mpfr_set_ui (z0, 0, MPFR_RNDN); 630 if (expected_inex > 0) 631 mpfr_nextabove (z0); 632 633 mpfr_clear_flags (); 634 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd); 635 cmpres (0, x, "emin - 2", (mpfr_rnd_t) rnd, z0, expected_inex, z, inex, 636 ufinex, "underflow_up3", "mpfr_pow"); 637 test_others (NULL, "emin - 2", (mpfr_rnd_t) rnd, x, y, z, inex, ufinex, 638 "underflow_up3"); 639 } 640 641 mpfr_clears (x, y, z, z0, (mpfr_ptr) 0); 642} 643 644static void 645underflow_up (void) 646{ 647 underflow_up1 (); 648 underflow_up2 (); 649 underflow_up3 (); 650} 651 652static void 653overflow_inv (void) 654{ 655 mpfr_t x, y, z; 656 int precx; 657 int s, t; 658 int inex; 659 int rnd; 660 661 mpfr_init2 (y, 2); 662 mpfr_init2 (z, 8); 663 664 mpfr_set_si (y, -1, MPFR_RNDN); 665 for (precx = 10; precx <= 100; precx += 90) 666 { 667 const char *sp = precx == 10 ? 668 "overflow_inv (precx = 10)" : "overflow_inv (precx = 100)"; 669 670 mpfr_init2 (x, precx); 671 for (s = -1; s <= 1; s += 2) 672 { 673 inex = mpfr_set_si_2exp (x, s, - mpfr_get_emax (), MPFR_RNDN); 674 MPFR_ASSERTN (inex == 0); 675 for (t = 0; t <= 5; t++) 676 { 677 /* If precx = 10: 678 * x = s * 2^(-emax) * (1 + t * 2^(-9)), so that 679 * 1/x = s * 2^emax * (1 - t * 2^(-9) + eps) with eps > 0. 680 * Values of (1/x) / 2^emax and overflow condition for x > 0: 681 * t = 0: 1 o: always 682 * t = 1: 0.11111111 100000000011... o: MPFR_RNDN and MPFR_RNDU 683 * t = 2: 0.11111111 000000001111... o: MPFR_RNDU 684 * t = 3: 0.11111110 100000100011... o: never 685 * 686 * If precx = 100: 687 * t = 0: always overflow 688 * t > 0: overflow for MPFR_RNDN and MPFR_RNDU. 689 */ 690 RND_LOOP_NO_RNDF (rnd) 691 { 692 int inf, overflow; 693 mpfr_rnd_t rnd2; 694 695 if (rnd == MPFR_RNDA) 696 rnd2 = s < 0 ? MPFR_RNDD : MPFR_RNDU; 697 else 698 rnd2 = (mpfr_rnd_t) rnd; 699 700 overflow = t == 0 || 701 ((mpfr_rnd_t) rnd == MPFR_RNDN && 702 (precx > 10 || t == 1)) || 703 (rnd2 == (s < 0 ? MPFR_RNDD : MPFR_RNDU) && 704 (precx > 10 || t <= 2)); 705 inf = overflow && 706 ((mpfr_rnd_t) rnd == MPFR_RNDN || 707 rnd2 == (s < 0 ? MPFR_RNDD : MPFR_RNDU)); 708 mpfr_clear_flags (); 709 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd); 710 if (overflow ^ !! mpfr_overflow_p ()) 711 { 712 printf ("Bad overflow flag in %s\nfor mpfr_pow%s\n" 713 "s = %d, t = %d, %s\n", sp, 714 ext ? ", extended exponent range" : "", 715 s, t, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 716 exit (1); 717 } 718 if (overflow && (inf ^ !! MPFR_IS_INF (z))) 719 { 720 printf ("Bad value in %s\nfor mpfr_pow%s\n" 721 "s = %d, t = %d, %s\nGot ", sp, 722 ext ? ", extended exponent range" : "", 723 s, t, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 724 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN); 725 printf (" instead of %s value.\n", 726 inf ? "infinite" : "finite"); 727 exit (1); 728 } 729 test_others (NULL, "-1", (mpfr_rnd_t) rnd, x, y, z, 730 inex, __gmpfr_flags, sp); 731 } /* RND_LOOP */ 732 mpfr_nexttoinf (x); 733 } /* for (t = ...) */ 734 } /* for (s = ...) */ 735 mpfr_clear (x); 736 } /* for (precx = ...) */ 737 738 mpfr_clears (y, z, (mpfr_ptr) 0); 739} 740 741static void 742alltst (void) 743{ 744 mpfr_exp_t emin, emax; 745 746 ext = 0; 747 tst (); 748 underflow_up (); 749 overflow_inv (); 750 751 emin = mpfr_get_emin (); 752 emax = mpfr_get_emax (); 753 set_emin (MPFR_EMIN_MIN); 754 set_emax (MPFR_EMAX_MAX); 755 if (mpfr_get_emin () != emin || mpfr_get_emax () != emax) 756 { 757 ext = 1; 758 tst (); 759 underflow_up (); 760 overflow_inv (); 761 set_emin (emin); 762 set_emax (emax); 763 } 764} 765 766int 767main (int argc, char *argv[]) 768{ 769 tests_start_mpfr (); 770 all_cmpres_errors = argc > 1; 771 alltst (); 772 tests_end_mpfr (); 773 return all_cmpres_errors < 0; 774} 775