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