1 /* Test file for mpfr_pow, mpfr_pow_ui and mpfr_pow_si. 2 3 Copyright 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 #include <stdio.h> 24 #include <stdlib.h> 25 #include <float.h> 26 #include <math.h> 27 #include <limits.h> 28 29 #include "mpfr-test.h" 30 31 #ifdef CHECK_EXTERNAL 32 static int 33 test_pow (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 34 { 35 int res; 36 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c) 37 && mpfr_get_prec (a) >= 53; 38 if (ok) 39 { 40 mpfr_print_raw (b); 41 printf (" "); 42 mpfr_print_raw (c); 43 } 44 res = mpfr_pow (a, b, c, rnd_mode); 45 if (ok) 46 { 47 printf (" "); 48 mpfr_print_raw (a); 49 printf ("\n"); 50 } 51 return res; 52 } 53 #else 54 #define test_pow mpfr_pow 55 #endif 56 57 #define TEST_FUNCTION test_pow 58 #define TWO_ARGS 59 #define TEST_RANDOM_POS 16 60 #define TGENERIC_NOWARNING 1 61 #include "tgeneric.c" 62 63 #define TEST_FUNCTION mpfr_pow_ui 64 #define INTEGER_TYPE unsigned long 65 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS) 66 #include "tgeneric_ui.c" 67 68 #define TEST_FUNCTION mpfr_pow_si 69 #define INTEGER_TYPE long 70 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS) 71 #define test_generic_ui test_generic_si 72 #include "tgeneric_ui.c" 73 74 static void 75 check_pow_ui (void) 76 { 77 mpfr_t a, b; 78 unsigned long n; 79 int res; 80 81 mpfr_init2 (a, 53); 82 mpfr_init2 (b, 53); 83 84 /* check in-place operations */ 85 mpfr_set_str (b, "0.6926773", 10, MPFR_RNDN); 86 mpfr_pow_ui (a, b, 10, MPFR_RNDN); 87 mpfr_pow_ui (b, b, 10, MPFR_RNDN); 88 if (mpfr_cmp (a, b)) 89 { 90 printf ("Error for mpfr_pow_ui (b, b, ...)\n"); 91 exit (1); 92 } 93 94 /* check large exponents */ 95 mpfr_set_ui (b, 1, MPFR_RNDN); 96 mpfr_pow_ui (a, b, 4294967295UL, MPFR_RNDN); 97 98 mpfr_set_inf (a, -1); 99 mpfr_pow_ui (a, a, 4049053855UL, MPFR_RNDN); 100 if (!mpfr_inf_p (a) || (mpfr_sgn (a) >= 0)) 101 { 102 printf ("Error for (-Inf)^4049053855\n"); 103 exit (1); 104 } 105 106 mpfr_set_inf (a, -1); 107 mpfr_pow_ui (a, a, (unsigned long) 30002752, MPFR_RNDN); 108 if (!mpfr_inf_p (a) || (mpfr_sgn (a) <= 0)) 109 { 110 printf ("Error for (-Inf)^30002752\n"); 111 exit (1); 112 } 113 114 /* Check underflow */ 115 mpfr_set_str_binary (a, "1E-1"); 116 res = mpfr_pow_ui (a, a, -mpfr_get_emin (), MPFR_RNDN); 117 if (MPFR_GET_EXP (a) != mpfr_get_emin () + 1) 118 { 119 printf ("Error for (1e-1)^MPFR_EMAX_MAX\n"); 120 mpfr_dump (a); 121 exit (1); 122 } 123 124 mpfr_set_str_binary (a, "1E-10"); 125 res = mpfr_pow_ui (a, a, -mpfr_get_emin (), MPFR_RNDZ); 126 if (!MPFR_IS_ZERO (a)) 127 { 128 printf ("Error for (1e-10)^MPFR_EMAX_MAX\n"); 129 mpfr_dump (a); 130 exit (1); 131 } 132 133 /* Check overflow */ 134 mpfr_set_str_binary (a, "1E10"); 135 res = mpfr_pow_ui (a, a, ULONG_MAX, MPFR_RNDN); 136 if (!MPFR_IS_INF (a) || MPFR_SIGN (a) < 0) 137 { 138 printf ("Error for (1e10)^ULONG_MAX\n"); 139 exit (1); 140 } 141 142 /* Bug in pow_ui.c from r3214 to r5107: if x = y (same mpfr_t argument), 143 the input argument is negative, n is odd, an overflow or underflow 144 occurs, and the temporary result res is positive, then the result 145 gets a wrong sign (positive instead of negative). */ 146 mpfr_set_str_binary (a, "-1E10"); 147 n = (ULONG_MAX ^ (ULONG_MAX >> 1)) + 1; 148 res = mpfr_pow_ui (a, a, n, MPFR_RNDN); 149 if (!MPFR_IS_INF (a) || MPFR_SIGN (a) > 0) 150 { 151 printf ("Error for (-1e10)^%lu, expected -Inf,\ngot ", n); 152 mpfr_dump (a); 153 exit (1); 154 } 155 156 /* Check 0 */ 157 MPFR_SET_ZERO (a); 158 MPFR_SET_POS (a); 159 mpfr_set_si (b, -1, MPFR_RNDN); 160 res = mpfr_pow_ui (b, a, 1, MPFR_RNDN); 161 if (res != 0 || MPFR_IS_NEG (b)) 162 { 163 printf ("Error for (0+)^1\n"); 164 exit (1); 165 } 166 MPFR_SET_ZERO (a); 167 MPFR_SET_NEG (a); 168 mpfr_set_ui (b, 1, MPFR_RNDN); 169 res = mpfr_pow_ui (b, a, 5, MPFR_RNDN); 170 if (res != 0 || MPFR_IS_POS (b)) 171 { 172 printf ("Error for (0-)^5\n"); 173 exit (1); 174 } 175 MPFR_SET_ZERO (a); 176 MPFR_SET_NEG (a); 177 mpfr_set_si (b, -1, MPFR_RNDN); 178 res = mpfr_pow_ui (b, a, 6, MPFR_RNDN); 179 if (res != 0 || MPFR_IS_NEG (b)) 180 { 181 printf ("Error for (0-)^6\n"); 182 exit (1); 183 } 184 185 mpfr_set_prec (a, 122); 186 mpfr_set_prec (b, 122); 187 mpfr_set_str_binary (a, "0.10000010010000111101001110100101101010011110011100001111000001001101000110011001001001001011001011010110110110101000111011E1"); 188 mpfr_set_str_binary (b, "0.11111111100101001001000001000001100011100000001110111111100011111000111011100111111111110100011000111011000100100011001011E51290375"); 189 mpfr_pow_ui (a, a, 2026876995UL, MPFR_RNDU); 190 if (mpfr_cmp (a, b) != 0) 191 { 192 printf ("Error for x^2026876995\n"); 193 exit (1); 194 } 195 196 mpfr_set_prec (a, 29); 197 mpfr_set_prec (b, 29); 198 mpfr_set_str_binary (a, "1.0000000000000000000000001111"); 199 mpfr_set_str_binary (b, "1.1001101111001100111001010111e165"); 200 mpfr_pow_ui (a, a, 2055225053, MPFR_RNDZ); 201 if (mpfr_cmp (a, b) != 0) 202 { 203 printf ("Error for x^2055225053\n"); 204 printf ("Expected "); 205 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN); 206 printf ("\nGot "); 207 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN); 208 printf ("\n"); 209 exit (1); 210 } 211 212 /* worst case found by Vincent Lefevre, 25 Nov 2006 */ 213 mpfr_set_prec (a, 53); 214 mpfr_set_prec (b, 53); 215 mpfr_set_str_binary (a, "1.0000010110000100001000101101101001011101101011010111"); 216 mpfr_set_str_binary (b, "1.0000110111101111011010110100001100010000001010110100E1"); 217 mpfr_pow_ui (a, a, 35, MPFR_RNDN); 218 if (mpfr_cmp (a, b) != 0) 219 { 220 printf ("Error in mpfr_pow_ui for worst case (1)\n"); 221 printf ("Expected "); 222 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN); 223 printf ("\nGot "); 224 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN); 225 printf ("\n"); 226 exit (1); 227 } 228 /* worst cases found on 2006-11-26 */ 229 mpfr_set_str_binary (a, "1.0110100111010001101001010111001110010100111111000011"); 230 mpfr_set_str_binary (b, "1.1111010011101110001111010110000101110000110110101100E17"); 231 mpfr_pow_ui (a, a, 36, MPFR_RNDD); 232 if (mpfr_cmp (a, b) != 0) 233 { 234 printf ("Error in mpfr_pow_ui for worst case (2)\n"); 235 printf ("Expected "); 236 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN); 237 printf ("\nGot "); 238 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN); 239 printf ("\n"); 240 exit (1); 241 } 242 mpfr_set_str_binary (a, "1.1001010100001110000110111111100011011101110011000100"); 243 mpfr_set_str_binary (b, "1.1100011101101101100010110001000001110001111110010001E23"); 244 mpfr_pow_ui (a, a, 36, MPFR_RNDU); 245 if (mpfr_cmp (a, b) != 0) 246 { 247 printf ("Error in mpfr_pow_ui for worst case (3)\n"); 248 printf ("Expected "); 249 mpfr_out_str (stdout, 2, 0, b, MPFR_RNDN); 250 printf ("\nGot "); 251 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN); 252 printf ("\n"); 253 exit (1); 254 } 255 256 mpfr_clear (a); 257 mpfr_clear (b); 258 } 259 260 static void 261 check_pow_si (void) 262 { 263 mpfr_t x; 264 265 mpfr_init (x); 266 267 mpfr_set_nan (x); 268 mpfr_pow_si (x, x, -1, MPFR_RNDN); 269 MPFR_ASSERTN(mpfr_nan_p (x)); 270 271 mpfr_set_inf (x, 1); 272 mpfr_pow_si (x, x, -1, MPFR_RNDN); 273 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x)); 274 275 mpfr_set_inf (x, -1); 276 mpfr_pow_si (x, x, -1, MPFR_RNDN); 277 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_NEG(x)); 278 279 mpfr_set_inf (x, -1); 280 mpfr_pow_si (x, x, -2, MPFR_RNDN); 281 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x)); 282 283 mpfr_set_ui (x, 0, MPFR_RNDN); 284 mpfr_pow_si (x, x, -1, MPFR_RNDN); 285 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0); 286 287 mpfr_set_ui (x, 0, MPFR_RNDN); 288 mpfr_neg (x, x, MPFR_RNDN); 289 mpfr_pow_si (x, x, -1, MPFR_RNDN); 290 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) < 0); 291 292 mpfr_set_ui (x, 0, MPFR_RNDN); 293 mpfr_neg (x, x, MPFR_RNDN); 294 mpfr_pow_si (x, x, -2, MPFR_RNDN); 295 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0); 296 297 mpfr_set_si (x, 2, MPFR_RNDN); 298 mpfr_pow_si (x, x, LONG_MAX, MPFR_RNDN); /* 2^LONG_MAX */ 299 if (LONG_MAX > mpfr_get_emax () - 1) /* LONG_MAX + 1 > emax */ 300 { 301 MPFR_ASSERTN (mpfr_inf_p (x)); 302 } 303 else 304 { 305 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, (mpfr_exp_t) LONG_MAX)); 306 } 307 308 mpfr_set_si (x, 2, MPFR_RNDN); 309 mpfr_pow_si (x, x, LONG_MIN, MPFR_RNDN); /* 2^LONG_MIN */ 310 if (LONG_MIN + 1 < mpfr_get_emin ()) 311 { 312 MPFR_ASSERTN (mpfr_zero_p (x)); 313 } 314 else 315 { 316 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, (mpfr_exp_t) LONG_MIN)); 317 } 318 319 mpfr_set_si (x, 2, MPFR_RNDN); 320 mpfr_pow_si (x, x, LONG_MIN + 1, MPFR_RNDN); /* 2^(LONG_MIN+1) */ 321 if (mpfr_nan_p (x)) 322 { 323 printf ("Error in pow_si(2, LONG_MIN+1): got NaN\n"); 324 exit (1); 325 } 326 if (LONG_MIN + 2 < mpfr_get_emin ()) 327 { 328 MPFR_ASSERTN (mpfr_zero_p (x)); 329 } 330 else 331 { 332 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, (mpfr_exp_t) (LONG_MIN + 1))); 333 } 334 335 mpfr_set_si_2exp (x, 1, -1, MPFR_RNDN); /* 0.5 */ 336 mpfr_pow_si (x, x, LONG_MIN, MPFR_RNDN); /* 2^(-LONG_MIN) */ 337 if (LONG_MIN < 1 - mpfr_get_emax ()) /* 1 - LONG_MIN > emax */ 338 { 339 MPFR_ASSERTN (mpfr_inf_p (x)); 340 } 341 else 342 { 343 MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 2, (mpfr_exp_t) - (LONG_MIN + 1))); 344 } 345 346 mpfr_clear (x); 347 } 348 349 static void 350 check_special_pow_si (void) 351 { 352 mpfr_t a, b; 353 mpfr_exp_t emin; 354 355 mpfr_init (a); 356 mpfr_init (b); 357 mpfr_set_str (a, "2E100000000", 10, MPFR_RNDN); 358 mpfr_set_si (b, -10, MPFR_RNDN); 359 test_pow (b, a, b, MPFR_RNDN); 360 if (!MPFR_IS_ZERO(b)) 361 { 362 printf("Pow(2E10000000, -10) failed\n"); 363 mpfr_dump (a); 364 mpfr_dump (b); 365 exit(1); 366 } 367 368 emin = mpfr_get_emin (); 369 mpfr_set_emin (-10); 370 mpfr_set_si (a, -2, MPFR_RNDN); 371 mpfr_pow_si (b, a, -10000, MPFR_RNDN); 372 if (!MPFR_IS_ZERO (b)) 373 { 374 printf ("Pow_so (1, -10000) doesn't underflow if emin=-10.\n"); 375 mpfr_dump (a); 376 mpfr_dump (b); 377 exit (1); 378 } 379 mpfr_set_emin (emin); 380 mpfr_clear (a); 381 mpfr_clear (b); 382 } 383 384 static void 385 pow_si_long_min (void) 386 { 387 mpfr_t x, y, z; 388 int inex; 389 390 mpfr_inits2 (sizeof(long) * CHAR_BIT + 32, x, y, z, (mpfr_ptr) 0); 391 mpfr_set_si_2exp (x, 3, -1, MPFR_RNDN); /* 1.5 */ 392 393 inex = mpfr_set_si (y, LONG_MIN, MPFR_RNDN); 394 MPFR_ASSERTN (inex == 0); 395 mpfr_nextbelow (y); 396 mpfr_pow (y, x, y, MPFR_RNDD); 397 398 inex = mpfr_set_si (z, LONG_MIN, MPFR_RNDN); 399 MPFR_ASSERTN (inex == 0); 400 mpfr_nextabove (z); 401 mpfr_pow (z, x, z, MPFR_RNDU); 402 403 mpfr_pow_si (x, x, LONG_MIN, MPFR_RNDN); /* 1.5^LONG_MIN */ 404 if (mpfr_cmp (x, y) < 0 || mpfr_cmp (x, z) > 0) 405 { 406 printf ("Error in pow_si_long_min\n"); 407 exit (1); 408 } 409 410 mpfr_clears (x, y, z, (mpfr_ptr) 0); 411 } 412 413 static void 414 check_inexact (mpfr_prec_t p) 415 { 416 mpfr_t x, y, z, t; 417 unsigned long u; 418 mpfr_prec_t q; 419 int inexact, cmp; 420 int rnd; 421 422 mpfr_init2 (x, p); 423 mpfr_init (y); 424 mpfr_init (z); 425 mpfr_init (t); 426 mpfr_urandomb (x, RANDS); 427 u = randlimb () % 2; 428 for (q = 2; q <= p; q++) 429 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 430 { 431 mpfr_set_prec (y, q); 432 mpfr_set_prec (z, q + 10); 433 mpfr_set_prec (t, q); 434 inexact = mpfr_pow_ui (y, x, u, (mpfr_rnd_t) rnd); 435 cmp = mpfr_pow_ui (z, x, u, (mpfr_rnd_t) rnd); 436 if (mpfr_can_round (z, q + 10, (mpfr_rnd_t) rnd, (mpfr_rnd_t) rnd, q)) 437 { 438 cmp = mpfr_set (t, z, (mpfr_rnd_t) rnd) || cmp; 439 if (mpfr_cmp (y, t)) 440 { 441 printf ("results differ for u=%lu rnd=%s\n", 442 u, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 443 printf ("x="); mpfr_print_binary (x); puts (""); 444 printf ("y="); mpfr_print_binary (y); puts (""); 445 printf ("t="); mpfr_print_binary (t); puts (""); 446 printf ("z="); mpfr_print_binary (z); puts (""); 447 exit (1); 448 } 449 if (((inexact == 0) && (cmp != 0)) || 450 ((inexact != 0) && (cmp == 0))) 451 { 452 printf ("Wrong inexact flag for p=%u, q=%u, rnd=%s\n", 453 (unsigned int) p, (unsigned int) q, 454 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 455 printf ("expected %d, got %d\n", cmp, inexact); 456 printf ("u=%lu x=", u); mpfr_print_binary (x); puts (""); 457 printf ("y="); mpfr_print_binary (y); puts (""); 458 exit (1); 459 } 460 } 461 } 462 463 /* check exact power */ 464 mpfr_set_prec (x, p); 465 mpfr_set_prec (y, p); 466 mpfr_set_prec (z, p); 467 mpfr_set_ui (x, 4, MPFR_RNDN); 468 mpfr_set_str (y, "0.5", 10, MPFR_RNDN); 469 test_pow (z, x, y, MPFR_RNDZ); 470 471 mpfr_clear (x); 472 mpfr_clear (y); 473 mpfr_clear (z); 474 mpfr_clear (t); 475 } 476 477 static void 478 special (void) 479 { 480 mpfr_t x, y, z, t; 481 mpfr_exp_t emin, emax; 482 int inex; 483 484 mpfr_init2 (x, 53); 485 mpfr_init2 (y, 53); 486 mpfr_init2 (z, 53); 487 mpfr_init2 (t, 2); 488 489 mpfr_set_ui (x, 2, MPFR_RNDN); 490 mpfr_pow_si (x, x, -2, MPFR_RNDN); 491 if (mpfr_cmp_ui_2exp (x, 1, -2)) 492 { 493 printf ("Error in pow_si(x,x,-2) for x=2\n"); 494 exit (1); 495 } 496 mpfr_set_ui (x, 2, MPFR_RNDN); 497 mpfr_set_si (y, -2, MPFR_RNDN); 498 test_pow (x, x, y, MPFR_RNDN); 499 if (mpfr_cmp_ui_2exp (x, 1, -2)) 500 { 501 printf ("Error in pow(x,x,y) for x=2, y=-2\n"); 502 exit (1); 503 } 504 505 mpfr_set_prec (x, 2); 506 mpfr_set_str_binary (x, "1.0e-1"); 507 mpfr_set_prec (y, 53); 508 mpfr_set_str_binary (y, "0.11010110011100101010110011001010100111000001000101110E-1"); 509 mpfr_set_prec (z, 2); 510 test_pow (z, x, y, MPFR_RNDZ); 511 mpfr_set_str_binary (x, "1.0e-1"); 512 if (mpfr_cmp (x, z)) 513 { 514 printf ("Error in mpfr_pow (1)\n"); 515 exit (1); 516 } 517 518 mpfr_set_prec (x, 64); 519 mpfr_set_prec (y, 64); 520 mpfr_set_prec (z, 64); 521 mpfr_set_prec (t, 64); 522 mpfr_set_str_binary (x, "0.111011000111100000111010000101010100110011010000011"); 523 mpfr_set_str_binary (y, "0.111110010100110000011101100011010111000010000100101"); 524 mpfr_set_str_binary (t, "0.1110110011110110001000110100100001001111010011111000010000011001"); 525 526 test_pow (z, x, y, MPFR_RNDN); 527 if (mpfr_cmp (z, t)) 528 { 529 printf ("Error in mpfr_pow for prec=64, rnd=MPFR_RNDN\n"); 530 exit (1); 531 } 532 533 mpfr_set_prec (x, 53); 534 mpfr_set_prec (y, 53); 535 mpfr_set_prec (z, 53); 536 mpfr_set_str (x, "5.68824667828621954868e-01", 10, MPFR_RNDN); 537 mpfr_set_str (y, "9.03327850535952658895e-01", 10, MPFR_RNDN); 538 test_pow (z, x, y, MPFR_RNDZ); 539 if (mpfr_cmp_str1 (z, "0.60071044650456473235")) 540 { 541 printf ("Error in mpfr_pow for prec=53, rnd=MPFR_RNDZ\n"); 542 exit (1); 543 } 544 545 mpfr_set_prec (t, 2); 546 mpfr_set_prec (x, 30); 547 mpfr_set_prec (y, 30); 548 mpfr_set_prec (z, 30); 549 mpfr_set_str (x, "1.00000000001010111110001111011e1", 2, MPFR_RNDN); 550 mpfr_set_str (t, "-0.5", 10, MPFR_RNDN); 551 test_pow (z, x, t, MPFR_RNDN); 552 mpfr_set_str (y, "1.01101001111010101110000101111e-1", 2, MPFR_RNDN); 553 if (mpfr_cmp (z, y)) 554 { 555 printf ("Error in mpfr_pow for prec=30, rnd=MPFR_RNDN\n"); 556 exit (1); 557 } 558 559 mpfr_set_prec (x, 21); 560 mpfr_set_prec (y, 21); 561 mpfr_set_prec (z, 21); 562 mpfr_set_str (x, "1.11111100100001100101", 2, MPFR_RNDN); 563 test_pow (z, x, t, MPFR_RNDZ); 564 mpfr_set_str (y, "1.01101011010001100000e-1", 2, MPFR_RNDN); 565 if (mpfr_cmp (z, y)) 566 { 567 printf ("Error in mpfr_pow for prec=21, rnd=MPFR_RNDZ\n"); 568 printf ("Expected "); 569 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 570 printf ("\nGot "); 571 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN); 572 printf ("\n"); 573 exit (1); 574 } 575 576 /* From http://www.terra.es/personal9/ismaeljc/hall.htm */ 577 mpfr_set_prec (x, 113); 578 mpfr_set_prec (y, 2); 579 mpfr_set_prec (z, 169); 580 mpfr_set_str1 (x, "6078673043126084065007902175846955"); 581 mpfr_set_ui_2exp (y, 3, -1, MPFR_RNDN); 582 test_pow (z, x, y, MPFR_RNDZ); 583 if (mpfr_cmp_str1 (z, "473928882491000966028828671876527456070714790264144")) 584 { 585 printf ("Error in mpfr_pow for 6078673043126084065007902175846955"); 586 printf ("^(3/2), MPFR_RNDZ\nExpected "); 587 printf ("4.73928882491000966028828671876527456070714790264144e50"); 588 printf ("\nGot "); 589 mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN); 590 printf ("\n"); 591 exit (1); 592 } 593 test_pow (z, x, y, MPFR_RNDU); 594 if (mpfr_cmp_str1 (z, "473928882491000966028828671876527456070714790264145")) 595 { 596 printf ("Error in mpfr_pow for 6078673043126084065007902175846955"); 597 printf ("^(3/2), MPFR_RNDU\nExpected "); 598 printf ("4.73928882491000966028828671876527456070714790264145e50"); 599 printf ("\nGot "); 600 mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN); 601 printf ("\n"); 602 exit (1); 603 } 604 605 mpfr_set_inf (x, 1); 606 mpfr_set_prec (y, 2); 607 mpfr_set_str_binary (y, "1E10"); 608 test_pow (z, x, y, MPFR_RNDN); 609 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z)); 610 mpfr_set_inf (x, -1); 611 test_pow (z, x, y, MPFR_RNDN); 612 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z)); 613 mpfr_set_prec (y, 10); 614 mpfr_set_str_binary (y, "1.000000001E9"); 615 test_pow (z, x, y, MPFR_RNDN); 616 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_NEG(z)); 617 mpfr_set_str_binary (y, "1.000000001E8"); 618 test_pow (z, x, y, MPFR_RNDN); 619 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z)); 620 621 mpfr_set_inf (x, -1); 622 mpfr_set_prec (y, 2 * mp_bits_per_limb); 623 mpfr_set_ui (y, 1, MPFR_RNDN); 624 mpfr_mul_2exp (y, y, mp_bits_per_limb - 1, MPFR_RNDN); 625 /* y = 2^(mp_bits_per_limb - 1) */ 626 test_pow (z, x, y, MPFR_RNDN); 627 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z)); 628 mpfr_nextabove (y); 629 test_pow (z, x, y, MPFR_RNDN); 630 /* y = 2^(mp_bits_per_limb - 1) + epsilon */ 631 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z)); 632 mpfr_nextbelow (y); 633 mpfr_div_2exp (y, y, 1, MPFR_RNDN); 634 mpfr_nextabove (y); 635 test_pow (z, x, y, MPFR_RNDN); 636 /* y = 2^(mp_bits_per_limb - 2) + epsilon */ 637 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z)); 638 639 mpfr_set_si (x, -1, MPFR_RNDN); 640 mpfr_set_prec (y, 2); 641 mpfr_set_str_binary (y, "1E10"); 642 test_pow (z, x, y, MPFR_RNDN); 643 MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0); 644 645 /* Check (-0)^(17.0001) */ 646 mpfr_set_prec (x, 6); 647 mpfr_set_prec (y, 640); 648 MPFR_SET_ZERO (x); MPFR_SET_NEG (x); 649 mpfr_set_ui (y, 17, MPFR_RNDN); mpfr_nextabove (y); 650 test_pow (z, x, y, MPFR_RNDN); 651 MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z)); 652 653 /* Bugs reported by Kevin Rauch on 29 Oct 2007 */ 654 emin = mpfr_get_emin (); 655 emax = mpfr_get_emax (); 656 mpfr_set_emin (-1000000); 657 mpfr_set_emax ( 1000000); 658 mpfr_set_prec (x, 64); 659 mpfr_set_prec (y, 64); 660 mpfr_set_prec (z, 64); 661 mpfr_set_str (x, "-0.5", 10, MPFR_RNDN); 662 mpfr_set_str (y, "-0.ffffffffffffffff", 16, MPFR_RNDN); 663 mpfr_set_exp (y, mpfr_get_emax ()); 664 inex = mpfr_pow (z, x, y, MPFR_RNDN); 665 /* (-0.5)^(-n) = 1/2^n for n even */ 666 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS (z) && inex > 0); 667 668 /* (-1)^(-n) = 1 for n even */ 669 mpfr_set_str (x, "-1", 10, MPFR_RNDN); 670 inex = mpfr_pow (z, x, y, MPFR_RNDN); 671 MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0 && inex == 0); 672 673 /* (-1)^n = 1 for n even */ 674 mpfr_set_str (x, "-1", 10, MPFR_RNDN); 675 mpfr_neg (y, y, MPFR_RNDN); 676 inex = mpfr_pow (z, x, y, MPFR_RNDN); 677 MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0 && inex == 0); 678 679 /* (-1.5)^n = +Inf for n even */ 680 mpfr_set_str (x, "-1.5", 10, MPFR_RNDN); 681 inex = mpfr_pow (z, x, y, MPFR_RNDN); 682 MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS (z) && inex > 0); 683 684 /* (-n)^1.5 = NaN for n even */ 685 mpfr_neg (y, y, MPFR_RNDN); 686 mpfr_set_str (x, "1.5", 10, MPFR_RNDN); 687 inex = mpfr_pow (z, y, x, MPFR_RNDN); 688 MPFR_ASSERTN(mpfr_nan_p (z)); 689 690 /* x^(-1.5) = NaN for x small < 0 */ 691 mpfr_set_str (x, "-0.8", 16, MPFR_RNDN); 692 mpfr_set_exp (x, mpfr_get_emin ()); 693 mpfr_set_str (y, "-1.5", 10, MPFR_RNDN); 694 inex = mpfr_pow (z, x, y, MPFR_RNDN); 695 MPFR_ASSERTN(mpfr_nan_p (z)); 696 697 mpfr_set_emin (emin); 698 mpfr_set_emax (emax); 699 mpfr_clear (x); 700 mpfr_clear (y); 701 mpfr_clear (z); 702 mpfr_clear (t); 703 } 704 705 static void 706 particular_cases (void) 707 { 708 mpfr_t t[11], r, r2; 709 mpz_t z; 710 long si; 711 712 static const char *name[11] = { 713 "NaN", "+inf", "-inf", "+0", "-0", "+1", "-1", "+2", "-2", "+0.5", "-0.5"}; 714 int i, j; 715 int error = 0; 716 717 mpz_init (z); 718 719 for (i = 0; i < 11; i++) 720 mpfr_init2 (t[i], 2); 721 mpfr_init2 (r, 6); 722 mpfr_init2 (r2, 6); 723 724 mpfr_set_nan (t[0]); 725 mpfr_set_inf (t[1], 1); 726 mpfr_set_ui (t[3], 0, MPFR_RNDN); 727 mpfr_set_ui (t[5], 1, MPFR_RNDN); 728 mpfr_set_ui (t[7], 2, MPFR_RNDN); 729 mpfr_div_2ui (t[9], t[5], 1, MPFR_RNDN); 730 for (i = 1; i < 11; i += 2) 731 mpfr_neg (t[i+1], t[i], MPFR_RNDN); 732 733 for (i = 0; i < 11; i++) 734 for (j = 0; j < 11; j++) 735 { 736 double d; 737 int p; 738 static const int q[11][11] = { 739 /* NaN +inf -inf +0 -0 +1 -1 +2 -2 +0.5 -0.5 */ 740 /* NaN */ { 0, 0, 0, 128, 128, 0, 0, 0, 0, 0, 0 }, 741 /* +inf */ { 0, 1, 2, 128, 128, 1, 2, 1, 2, 1, 2 }, 742 /* -inf */ { 0, 1, 2, 128, 128, -1, -2, 1, 2, 1, 2 }, 743 /* +0 */ { 0, 2, 1, 128, 128, 2, 1, 2, 1, 2, 1 }, 744 /* -0 */ { 0, 2, 1, 128, 128, -2, -1, 2, 1, 2, 1 }, 745 /* +1 */ {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128 }, 746 /* -1 */ { 0, 128, 128, 128, 128,-128,-128, 128, 128, 0, 0 }, 747 /* +2 */ { 0, 1, 2, 128, 128, 256, 64, 512, 32, 180, 90 }, 748 /* -2 */ { 0, 1, 2, 128, 128,-256, -64, 512, 32, 0, 0 }, 749 /* +0.5 */ { 0, 2, 1, 128, 128, 64, 256, 32, 512, 90, 180 }, 750 /* -0.5 */ { 0, 2, 1, 128, 128, -64,-256, 32, 512, 0, 0 } 751 }; 752 /* This define is used to make the following table readable */ 753 #define N MPFR_FLAGS_NAN 754 #define I MPFR_FLAGS_INEXACT 755 #define D MPFR_FLAGS_DIVBY0 756 static const unsigned int f[11][11] = { 757 /* NaN +inf -inf +0 -0 +1 -1 +2 -2 +0.5 -0.5 */ 758 /* NaN */ { N, N, N, 0, 0, N, N, N, N, N, N }, 759 /* +inf */ { N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, 760 /* -inf */ { N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, 761 /* +0 */ { N, 0, 0, 0, 0, 0, D, 0, D, 0, D }, 762 /* -0 */ { N, 0, 0, 0, 0, 0, D, 0, D, 0, D }, 763 /* +1 */ { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, 764 /* -1 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, N, N }, 765 /* +2 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, I, I }, 766 /* -2 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, N, N }, 767 /* +0.5 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, I, I }, 768 /* -0.5 */ { N, 0, 0, 0, 0, 0, 0, 0, 0, N, N } 769 }; 770 #undef N 771 #undef I 772 #undef D 773 mpfr_clear_flags (); 774 test_pow (r, t[i], t[j], MPFR_RNDN); 775 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 : 776 mpfr_cmp_ui (r, 0) == 0 ? 2 : 777 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0)); 778 if (p != 0 && MPFR_IS_NEG (r)) 779 p = -p; 780 if (p != q[i][j]) 781 { 782 printf ("Error in mpfr_pow for (%s)^(%s) (%d,%d):\n" 783 "got %d instead of %d\n", 784 name[i], name[j], i, j, p, q[i][j]); 785 mpfr_dump (r); 786 error = 1; 787 } 788 if (__gmpfr_flags != f[i][j]) 789 { 790 printf ("Error in mpfr_pow for (%s)^(%s) (%d,%d):\n" 791 "Flags = %u instead of expected %u\n", 792 name[i], name[j], i, j, __gmpfr_flags, f[i][j]); 793 mpfr_dump (r); 794 error = 1; 795 } 796 /* Perform the same tests with pow_z & pow_si & pow_ui 797 if t[j] is an integer */ 798 if (mpfr_integer_p (t[j])) 799 { 800 /* mpfr_pow_z */ 801 mpfr_clear_flags (); 802 mpfr_get_z (z, t[j], MPFR_RNDN); 803 mpfr_pow_z (r, t[i], z, MPFR_RNDN); 804 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 : 805 mpfr_cmp_ui (r, 0) == 0 ? 2 : 806 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0)); 807 if (p != 0 && MPFR_IS_NEG (r)) 808 p = -p; 809 if (p != q[i][j]) 810 { 811 printf ("Error in mpfr_pow_z for (%s)^(%s) (%d,%d):\n" 812 "got %d instead of %d\n", 813 name[i], name[j], i, j, p, q[i][j]); 814 mpfr_dump (r); 815 error = 1; 816 } 817 if (__gmpfr_flags != f[i][j]) 818 { 819 printf ("Error in mpfr_pow_z for (%s)^(%s) (%d,%d):\n" 820 "Flags = %u instead of expected %u\n", 821 name[i], name[j], i, j, __gmpfr_flags, f[i][j]); 822 mpfr_dump (r); 823 error = 1; 824 } 825 /* mpfr_pow_si */ 826 mpfr_clear_flags (); 827 si = mpfr_get_si (t[j], MPFR_RNDN); 828 mpfr_pow_si (r, t[i], si, MPFR_RNDN); 829 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 : 830 mpfr_cmp_ui (r, 0) == 0 ? 2 : 831 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0)); 832 if (p != 0 && MPFR_IS_NEG (r)) 833 p = -p; 834 if (p != q[i][j]) 835 { 836 printf ("Error in mpfr_pow_si for (%s)^(%s) (%d,%d):\n" 837 "got %d instead of %d\n", 838 name[i], name[j], i, j, p, q[i][j]); 839 mpfr_dump (r); 840 error = 1; 841 } 842 if (__gmpfr_flags != f[i][j]) 843 { 844 printf ("Error in mpfr_pow_si for (%s)^(%s) (%d,%d):\n" 845 "Flags = %u instead of expected %u\n", 846 name[i], name[j], i, j, __gmpfr_flags, f[i][j]); 847 mpfr_dump (r); 848 error = 1; 849 } 850 /* if si >= 0, test mpfr_pow_ui */ 851 if (si >= 0) 852 { 853 mpfr_clear_flags (); 854 mpfr_pow_ui (r, t[i], si, MPFR_RNDN); 855 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 : 856 mpfr_cmp_ui (r, 0) == 0 ? 2 : 857 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0)); 858 if (p != 0 && MPFR_IS_NEG (r)) 859 p = -p; 860 if (p != q[i][j]) 861 { 862 printf ("Error in mpfr_pow_ui for (%s)^(%s) (%d,%d):\n" 863 "got %d instead of %d\n", 864 name[i], name[j], i, j, p, q[i][j]); 865 mpfr_dump (r); 866 error = 1; 867 } 868 if (__gmpfr_flags != f[i][j]) 869 { 870 printf ("Error in mpfr_pow_ui for (%s)^(%s) (%d,%d):\n" 871 "Flags = %u instead of expected %u\n", 872 name[i], name[j], i, j, __gmpfr_flags, f[i][j]); 873 mpfr_dump (r); 874 error = 1; 875 } 876 } 877 } /* integer_p */ 878 /* Perform the same tests with mpfr_ui_pow */ 879 if (mpfr_integer_p (t[i]) && MPFR_IS_POS (t[i])) 880 { 881 /* mpfr_ui_pow */ 882 mpfr_clear_flags (); 883 si = mpfr_get_si (t[i], MPFR_RNDN); 884 mpfr_ui_pow (r, si, t[j], MPFR_RNDN); 885 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 : 886 mpfr_cmp_ui (r, 0) == 0 ? 2 : 887 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0)); 888 if (p != 0 && MPFR_IS_NEG (r)) 889 p = -p; 890 if (p != q[i][j]) 891 { 892 printf ("Error in mpfr_ui_pow for (%s)^(%s) (%d,%d):\n" 893 "got %d instead of %d\n", 894 name[i], name[j], i, j, p, q[i][j]); 895 mpfr_dump (r); 896 error = 1; 897 } 898 if (__gmpfr_flags != f[i][j]) 899 { 900 printf ("Error in mpfr_ui_pow for (%s)^(%s) (%d,%d):\n" 901 "Flags = %u instead of expected %u\n", 902 name[i], name[j], i, j, __gmpfr_flags, f[i][j]); 903 mpfr_dump (r); 904 error = 1; 905 } 906 } 907 } 908 909 for (i = 0; i < 11; i++) 910 mpfr_clear (t[i]); 911 mpfr_clear (r); 912 mpfr_clear (r2); 913 mpz_clear (z); 914 915 if (error) 916 exit (1); 917 } 918 919 static void 920 underflows (void) 921 { 922 mpfr_t x, y, z; 923 int err = 0; 924 int inexact; 925 int i; 926 mpfr_exp_t emin; 927 928 mpfr_init2 (x, 64); 929 mpfr_init2 (y, 64); 930 931 mpfr_set_ui (x, 1, MPFR_RNDN); 932 mpfr_set_exp (x, mpfr_get_emin()); 933 934 for (i = 3; i < 10; i++) 935 { 936 mpfr_set_ui (y, i, MPFR_RNDN); 937 mpfr_div_2ui (y, y, 1, MPFR_RNDN); 938 test_pow (y, x, y, MPFR_RNDN); 939 if (!MPFR_IS_FP(y) || mpfr_cmp_ui (y, 0)) 940 { 941 printf ("Error in mpfr_pow for "); 942 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 943 printf (" ^ (%d/2)\nGot ", i); 944 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 945 printf (" instead of 0.\n"); 946 exit (1); 947 } 948 } 949 950 mpfr_init2 (z, 55); 951 mpfr_set_str (x, "0.110011010011101001110001110100010000110111101E0", 952 2, MPFR_RNDN); 953 mpfr_set_str (y, "0.101110010011111001011010100011011100111110011E40", 954 2, MPFR_RNDN); 955 mpfr_clear_flags (); 956 inexact = mpfr_pow (z, x, y, MPFR_RNDU); 957 if (!mpfr_underflow_p ()) 958 { 959 printf ("Underflow flag is not set for special underflow test.\n"); 960 err = 1; 961 } 962 if (inexact <= 0) 963 { 964 printf ("Ternary value is wrong for special underflow test.\n"); 965 err = 1; 966 } 967 mpfr_set_ui (x, 0, MPFR_RNDN); 968 mpfr_nextabove (x); 969 if (mpfr_cmp (x, z) != 0) 970 { 971 printf ("Wrong value for special underflow test.\nGot "); 972 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN); 973 printf ("\ninstead of "); 974 mpfr_out_str (stdout, 2, 2, x, MPFR_RNDN); 975 printf ("\n"); 976 err = 1; 977 } 978 if (err) 979 exit (1); 980 981 /* MPFR currently (2006-08-19) segfaults on the following code (and 982 possibly makes other programs crash due to the lack of memory), 983 because y is converted into an mpz_t, and the required precision 984 is too high. */ 985 mpfr_set_prec (x, 2); 986 mpfr_set_prec (y, 2); 987 mpfr_set_prec (z, 12); 988 mpfr_set_ui_2exp (x, 3, -2, MPFR_RNDN); 989 mpfr_set_ui_2exp (y, 1, mpfr_get_emax () - 1, MPFR_RNDN); 990 mpfr_clear_flags (); 991 mpfr_pow (z, x, y, MPFR_RNDN); 992 if (!mpfr_underflow_p () || MPFR_NOTZERO (z)) 993 { 994 printf ("Underflow test with large y fails.\n"); 995 exit (1); 996 } 997 998 emin = mpfr_get_emin (); 999 mpfr_set_emin (-256); 1000 mpfr_set_prec (x, 2); 1001 mpfr_set_prec (y, 2); 1002 mpfr_set_prec (z, 12); 1003 mpfr_set_ui_2exp (x, 3, -2, MPFR_RNDN); 1004 mpfr_set_ui_2exp (y, 1, 38, MPFR_RNDN); 1005 mpfr_clear_flags (); 1006 inexact = mpfr_pow (z, x, y, MPFR_RNDN); 1007 if (!mpfr_underflow_p () || MPFR_NOTZERO (z) || inexact >= 0) 1008 { 1009 printf ("Bad underflow detection for 0.75^(2^38). Obtained:\n" 1010 "Underflow flag... %-3s (should be 'yes')\n" 1011 "Zero result...... %-3s (should be 'yes')\n" 1012 "Inexact value.... %-3d (should be negative)\n", 1013 mpfr_underflow_p () ? "yes" : "no", 1014 MPFR_IS_ZERO (z) ? "yes" : "no", inexact); 1015 exit (1); 1016 } 1017 mpfr_set_emin (emin); 1018 1019 emin = mpfr_get_emin (); 1020 mpfr_set_emin (-256); 1021 mpfr_set_prec (x, 2); 1022 mpfr_set_prec (y, 40); 1023 mpfr_set_prec (z, 12); 1024 mpfr_set_ui_2exp (x, 3, -1, MPFR_RNDN); 1025 mpfr_set_si_2exp (y, -1, 38, MPFR_RNDN); 1026 for (i = 0; i < 4; i++) 1027 { 1028 if (i == 2) 1029 mpfr_neg (x, x, MPFR_RNDN); 1030 mpfr_clear_flags (); 1031 inexact = mpfr_pow (z, x, y, MPFR_RNDN); 1032 if (!mpfr_underflow_p () || MPFR_NOTZERO (z) || 1033 (i == 3 ? (inexact <= 0) : (inexact >= 0))) 1034 { 1035 printf ("Bad underflow detection for ("); 1036 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 1037 printf (")^(-2^38-%d). Obtained:\n" 1038 "Overflow flag.... %-3s (should be 'no')\n" 1039 "Underflow flag... %-3s (should be 'yes')\n" 1040 "Zero result...... %-3s (should be 'yes')\n" 1041 "Inexact value.... %-3d (should be %s)\n", i, 1042 mpfr_overflow_p () ? "yes" : "no", 1043 mpfr_underflow_p () ? "yes" : "no", 1044 MPFR_IS_ZERO (z) ? "yes" : "no", inexact, 1045 i == 3 ? "positive" : "negative"); 1046 exit (1); 1047 } 1048 inexact = mpfr_sub_ui (y, y, 1, MPFR_RNDN); 1049 MPFR_ASSERTN (inexact == 0); 1050 } 1051 mpfr_set_emin (emin); 1052 1053 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1054 } 1055 1056 static void 1057 overflows (void) 1058 { 1059 mpfr_t a, b; 1060 1061 /* bug found by Ming J. Tsai <mingjt@delvron.us>, 4 Oct 2003 */ 1062 1063 mpfr_init_set_str (a, "5.1e32", 10, MPFR_RNDN); 1064 mpfr_init (b); 1065 1066 test_pow (b, a, a, MPFR_RNDN); 1067 if (!(mpfr_inf_p (b) && mpfr_sgn (b) > 0)) 1068 { 1069 printf ("Error for a^a for a=5.1e32\n"); 1070 printf ("Expected +Inf, got "); 1071 mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN); 1072 printf ("\n"); 1073 exit (1); 1074 } 1075 1076 mpfr_clear(a); 1077 mpfr_clear(b); 1078 } 1079 1080 static void 1081 overflows2 (void) 1082 { 1083 mpfr_t x, y, z; 1084 mpfr_exp_t emin, emax; 1085 int e; 1086 1087 /* x^y in reduced exponent range, where x = 2^b and y is not an integer 1088 (so that mpfr_pow_z is not used). */ 1089 1090 emin = mpfr_get_emin (); 1091 emax = mpfr_get_emax (); 1092 set_emin (-128); 1093 1094 mpfr_inits2 (16, x, y, z, (mpfr_ptr) 0); 1095 1096 mpfr_set_si_2exp (x, 1, -64, MPFR_RNDN); /* 2^(-64) */ 1097 mpfr_set_si_2exp (y, -1, -1, MPFR_RNDN); /* -0.5 */ 1098 for (e = 2; e <= 32; e += 17) 1099 { 1100 set_emax (e); 1101 mpfr_clear_flags (); 1102 mpfr_pow (z, x, y, MPFR_RNDN); 1103 if (MPFR_IS_NEG (z) || ! mpfr_inf_p (z)) 1104 { 1105 printf ("Error in overflows2 (e = %d): expected +Inf, got ", e); 1106 mpfr_dump (z); 1107 exit (1); 1108 } 1109 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 1110 { 1111 printf ("Error in overflows2 (e = %d): bad flags (%u)\n", 1112 e, __gmpfr_flags); 1113 exit (1); 1114 } 1115 } 1116 1117 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1118 1119 set_emin (emin); 1120 set_emax (emax); 1121 } 1122 1123 static void 1124 overflows3 (void) 1125 { 1126 /* x^y where x = 2^b, y is not an integer (so that mpfr_pow_z is not used) 1127 and b * y = emax in the extended exponent range. If emax is divisible 1128 by 3, we choose x = 2^(-2*emax/3) and y = -3/2. 1129 Test also with nextbelow(x). */ 1130 1131 if (MPFR_EMAX_MAX % 3 == 0) 1132 { 1133 mpfr_t x, y, z, t; 1134 mpfr_exp_t emin, emax; 1135 unsigned int flags; 1136 int i; 1137 1138 emin = mpfr_get_emin (); 1139 emax = mpfr_get_emax (); 1140 set_emin (MPFR_EMIN_MIN); 1141 set_emax (MPFR_EMAX_MAX); 1142 1143 mpfr_inits2 (16, x, y, z, t, (mpfr_ptr) 0); 1144 1145 mpfr_set_si_2exp (x, 1, -2 * (MPFR_EMAX_MAX / 3), MPFR_RNDN); 1146 for (i = 0; i <= 1; i++) 1147 { 1148 mpfr_set_si_2exp (y, -3, -1, MPFR_RNDN); 1149 mpfr_clear_flags (); 1150 mpfr_pow (z, x, y, MPFR_RNDN); 1151 if (MPFR_IS_NEG (z) || ! mpfr_inf_p (z)) 1152 { 1153 printf ("Error in overflows3 (RNDN, i = %d): expected +Inf," 1154 " got ", i); 1155 mpfr_dump (z); 1156 exit (1); 1157 } 1158 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 1159 { 1160 printf ("Error in overflows3 (RNDN, i = %d): bad flags (%u)\n", 1161 i, __gmpfr_flags); 1162 exit (1); 1163 } 1164 1165 mpfr_clear_flags (); 1166 mpfr_pow (z, x, y, MPFR_RNDZ); 1167 flags = __gmpfr_flags; 1168 mpfr_set (t, z, MPFR_RNDN); 1169 mpfr_nextabove (t); 1170 if (MPFR_IS_NEG (z) || mpfr_inf_p (z) || ! mpfr_inf_p (t)) 1171 { 1172 printf ("Error in overflows3 (RNDZ, i = %d):\nexpected ", i); 1173 mpfr_set_inf (t, 1); 1174 mpfr_nextbelow (t); 1175 mpfr_dump (t); 1176 printf ("got "); 1177 mpfr_dump (z); 1178 exit (1); 1179 } 1180 if (flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 1181 { 1182 printf ("Error in overflows3 (RNDZ, i = %d): bad flags (%u)\n", 1183 i, flags); 1184 exit (1); 1185 } 1186 mpfr_nextbelow (x); 1187 } 1188 1189 mpfr_clears (x, y, z, t, (mpfr_ptr) 0); 1190 1191 set_emin (emin); 1192 set_emax (emax); 1193 } 1194 } 1195 1196 static void 1197 x_near_one (void) 1198 { 1199 mpfr_t x, y, z; 1200 int inex; 1201 1202 mpfr_init2 (x, 32); 1203 mpfr_init2 (y, 4); 1204 mpfr_init2 (z, 33); 1205 1206 mpfr_set_ui (x, 1, MPFR_RNDN); 1207 mpfr_nextbelow (x); 1208 mpfr_set_ui_2exp (y, 11, -2, MPFR_RNDN); 1209 inex = mpfr_pow (z, x, y, MPFR_RNDN); 1210 if (mpfr_cmp_str (z, "0.111111111111111111111111111111011E0", 2, MPFR_RNDN) 1211 || inex <= 0) 1212 { 1213 printf ("Failure in x_near_one, got inex = %d and\nz = ", inex); 1214 mpfr_dump (z); 1215 } 1216 1217 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1218 } 1219 1220 static int 1221 mpfr_pow275 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) 1222 { 1223 mpfr_t z; 1224 int inex; 1225 1226 mpfr_init2 (z, 4); 1227 mpfr_set_ui_2exp (z, 11, -2, MPFR_RNDN); 1228 inex = mpfr_pow (y, x, z, MPFR_RNDN); 1229 mpfr_clear (z); 1230 return inex; 1231 } 1232 1233 /* Bug found by Kevin P. Rauch */ 1234 static void 1235 bug20071103 (void) 1236 { 1237 mpfr_t x, y, z; 1238 mpfr_exp_t emin, emax; 1239 1240 emin = mpfr_get_emin (); 1241 emax = mpfr_get_emax (); 1242 mpfr_set_emin (-1000000); 1243 mpfr_set_emax ( 1000000); 1244 1245 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0); 1246 mpfr_set_si_2exp (x, -3, -1, MPFR_RNDN); /* x = -1.5 */ 1247 mpfr_set_str (y, "-0.ffffffffffffffff", 16, MPFR_RNDN); 1248 mpfr_set_exp (y, mpfr_get_emax ()); 1249 mpfr_clear_flags (); 1250 mpfr_pow (z, x, y, MPFR_RNDN); 1251 MPFR_ASSERTN (mpfr_zero_p (z) && MPFR_SIGN (z) > 0 && 1252 __gmpfr_flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT)); 1253 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1254 1255 set_emin (emin); 1256 set_emax (emax); 1257 } 1258 1259 /* Bug found by Kevin P. Rauch */ 1260 static void 1261 bug20071104 (void) 1262 { 1263 mpfr_t x, y, z; 1264 mpfr_exp_t emin, emax; 1265 int inex; 1266 1267 emin = mpfr_get_emin (); 1268 emax = mpfr_get_emax (); 1269 mpfr_set_emin (-1000000); 1270 mpfr_set_emax ( 1000000); 1271 1272 mpfr_inits2 (20, x, y, z, (mpfr_ptr) 0); 1273 mpfr_set_ui (x, 0, MPFR_RNDN); 1274 mpfr_nextbelow (x); /* x = -2^(emin-1) */ 1275 mpfr_set_si (y, -2, MPFR_RNDN); /* y = -2 */ 1276 mpfr_clear_flags (); 1277 inex = mpfr_pow (z, x, y, MPFR_RNDN); 1278 if (! mpfr_inf_p (z) || MPFR_SIGN (z) < 0) 1279 { 1280 printf ("Error in bug20071104: expected +Inf, got "); 1281 mpfr_dump (z); 1282 exit (1); 1283 } 1284 if (inex <= 0) 1285 { 1286 printf ("Error in bug20071104: bad ternary value (%d)\n", inex); 1287 exit (1); 1288 } 1289 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 1290 { 1291 printf ("Error in bug20071104: bad flags (%u)\n", __gmpfr_flags); 1292 exit (1); 1293 } 1294 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1295 1296 set_emin (emin); 1297 set_emax (emax); 1298 } 1299 1300 /* Bug found by Kevin P. Rauch */ 1301 static void 1302 bug20071127 (void) 1303 { 1304 mpfr_t x, y, z; 1305 int i, tern; 1306 mpfr_exp_t emin, emax; 1307 1308 emin = mpfr_get_emin (); 1309 emax = mpfr_get_emax (); 1310 mpfr_set_emin (-1000000); 1311 mpfr_set_emax ( 1000000); 1312 1313 mpfr_init2 (x, 128); 1314 mpfr_init2 (y, 128); 1315 mpfr_init2 (z, 128); 1316 1317 mpfr_set_str (x, "0.80000000000000000000000000000001", 16, MPFR_RNDN); 1318 1319 for (i = 1; i < 9; i *= 2) 1320 { 1321 mpfr_set_str (y, "8000000000000000", 16, MPFR_RNDN); 1322 mpfr_add_si (y, y, i, MPFR_RNDN); 1323 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1324 MPFR_ASSERTN (mpfr_zero_p (z) && MPFR_IS_POS (z) && tern < 0); 1325 } 1326 1327 mpfr_clear (x); 1328 mpfr_clear (y); 1329 mpfr_clear (z); 1330 1331 mpfr_set_emin (emin); 1332 mpfr_set_emax (emax); 1333 } 1334 1335 /* Bug found by Kevin P. Rauch */ 1336 static void 1337 bug20071128 (void) 1338 { 1339 mpfr_t max_val, x, y, z; 1340 int i, tern; 1341 mpfr_exp_t emin, emax; 1342 1343 emin = mpfr_get_emin (); 1344 emax = mpfr_get_emax (); 1345 mpfr_set_emin (-1000000); 1346 mpfr_set_emax ( 1000000); 1347 1348 mpfr_init2 (max_val, 64); 1349 mpfr_init2 (x, 64); 1350 mpfr_init2 (y, 64); 1351 mpfr_init2 (z, 64); 1352 1353 mpfr_set_str (max_val, "0.ffffffffffffffff", 16, MPFR_RNDN); 1354 mpfr_set_exp (max_val, mpfr_get_emax ()); 1355 1356 mpfr_neg (x, max_val, MPFR_RNDN); 1357 1358 /* on 64-bit machines */ 1359 for (i = 41; i < 45; i++) 1360 { 1361 mpfr_set_si_2exp (y, -1, i, MPFR_RNDN); 1362 mpfr_add_si (y, y, 1, MPFR_RNDN); 1363 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1364 MPFR_ASSERTN (mpfr_zero_p (z) && MPFR_IS_NEG (z) && tern > 0); 1365 } 1366 1367 /* on 32-bit machines */ 1368 for (i = 9; i < 13; i++) 1369 { 1370 mpfr_set_si_2exp (y, -1, i, MPFR_RNDN); 1371 mpfr_add_si (y, y, 1, MPFR_RNDN); 1372 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1373 MPFR_ASSERTN(mpfr_zero_p (z) && MPFR_SIGN(z) < 0); 1374 } 1375 1376 mpfr_clear (x); 1377 mpfr_clear (y); 1378 mpfr_clear (z); 1379 mpfr_clear (max_val); 1380 1381 mpfr_set_emin (emin); 1382 mpfr_set_emax (emax); 1383 } 1384 1385 /* Bug found by Kevin P. Rauch */ 1386 static void 1387 bug20071218 (void) 1388 { 1389 mpfr_t x, y, z, t; 1390 int tern; 1391 1392 mpfr_inits2 (64, x, y, z, t, (mpfr_ptr) 0); 1393 mpfr_set_str (x, "0x.80000000000002P-1023", 0, MPFR_RNDN); 1394 mpfr_set_str (y, "100000.000000002", 16, MPFR_RNDN); 1395 mpfr_set_ui (t, 0, MPFR_RNDN); 1396 mpfr_nextabove (t); 1397 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1398 if (mpfr_cmp0 (z, t) != 0) 1399 { 1400 printf ("Error in bug20071218 (1): Expected\n"); 1401 mpfr_dump (t); 1402 printf ("Got\n"); 1403 mpfr_dump (z); 1404 exit (1); 1405 } 1406 if (tern <= 0) 1407 { 1408 printf ("Error in bug20071218 (1): bad ternary value" 1409 " (%d instead of positive)\n", tern); 1410 exit (1); 1411 } 1412 mpfr_mul_2ui (y, y, 32, MPFR_RNDN); 1413 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1414 if (MPFR_NOTZERO (z) || MPFR_IS_NEG (z)) 1415 { 1416 printf ("Error in bug20071218 (2): expected 0, got\n"); 1417 mpfr_dump (z); 1418 exit (1); 1419 } 1420 if (tern >= 0) 1421 { 1422 printf ("Error in bug20071218 (2): bad ternary value" 1423 " (%d instead of negative)\n", tern); 1424 exit (1); 1425 } 1426 mpfr_clears (x, y, z, t, (mpfr_ptr) 0); 1427 } 1428 1429 /* With revision 5429, this gives: 1430 * pow.c:43: assertion failed: !mpfr_integer_p (y) 1431 * This is fixed in revision 5432. 1432 */ 1433 static void 1434 bug20080721 (void) 1435 { 1436 mpfr_t x, y, z, t[2]; 1437 int inex; 1438 int rnd; 1439 int err = 0; 1440 1441 /* Note: input values have been chosen in a way to select the 1442 * general case. If mpfr_pow is modified, in particular line 1443 * if (y_is_integer && (MPFR_GET_EXP (y) <= 256)) 1444 * make sure that this test still does what we want. 1445 */ 1446 mpfr_inits2 (4913, x, y, (mpfr_ptr) 0); 1447 mpfr_inits2 (8, z, t[0], t[1], (mpfr_ptr) 0); 1448 mpfr_set_si (x, -1, MPFR_RNDN); 1449 mpfr_nextbelow (x); 1450 mpfr_set_ui_2exp (y, 1, mpfr_get_prec (y) - 1, MPFR_RNDN); 1451 inex = mpfr_add_ui (y, y, 1, MPFR_RNDN); 1452 MPFR_ASSERTN (inex == 0); 1453 mpfr_set_str_binary (t[0], "-0.10101101e2"); 1454 mpfr_set_str_binary (t[1], "-0.10101110e2"); 1455 RND_LOOP (rnd) 1456 { 1457 int i, inex0; 1458 1459 i = (rnd == MPFR_RNDN || rnd == MPFR_RNDD || rnd == MPFR_RNDA); 1460 inex0 = i ? -1 : 1; 1461 mpfr_clear_flags (); 1462 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd); 1463 if (__gmpfr_flags != MPFR_FLAGS_INEXACT || ! SAME_SIGN (inex, inex0) 1464 || MPFR_IS_NAN (z) || mpfr_cmp (z, t[i]) != 0) 1465 { 1466 unsigned int flags = __gmpfr_flags; 1467 1468 printf ("Error in bug20080721 with %s\n", 1469 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 1470 printf ("expected "); 1471 mpfr_out_str (stdout, 2, 0, t[i], MPFR_RNDN); 1472 printf (", inex = %d, flags = %u\n", inex0, 1473 (unsigned int) MPFR_FLAGS_INEXACT); 1474 printf ("got "); 1475 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN); 1476 printf (", inex = %d, flags = %u\n", inex, flags); 1477 err = 1; 1478 } 1479 } 1480 mpfr_clears (x, y, z, t[0], t[1], (mpfr_ptr) 0); 1481 if (err) 1482 exit (1); 1483 } 1484 1485 /* The following test fails in r5552 (32-bit and 64-bit). This is due to: 1486 * mpfr_log (t, absx, MPFR_RNDU); 1487 * mpfr_mul (t, y, t, MPFR_RNDU); 1488 * in pow.c, that is supposed to compute an upper bound on exp(y*ln|x|), 1489 * but this is incorrect if y is negative. 1490 */ 1491 static void 1492 bug20080820 (void) 1493 { 1494 mpfr_exp_t emin; 1495 mpfr_t x, y, z1, z2; 1496 1497 emin = mpfr_get_emin (); 1498 mpfr_set_emin (MPFR_EMIN_MIN); 1499 mpfr_init2 (x, 80); 1500 mpfr_init2 (y, sizeof (mpfr_exp_t) * CHAR_BIT + 32); 1501 mpfr_init2 (z1, 2); 1502 mpfr_init2 (z2, 80); 1503 mpfr_set_ui (x, 2, MPFR_RNDN); 1504 mpfr_nextbelow (x); 1505 mpfr_set_exp_t (y, mpfr_get_emin () - 2, MPFR_RNDN); 1506 mpfr_nextabove (y); 1507 mpfr_pow (z1, x, y, MPFR_RNDN); 1508 mpfr_pow (z2, x, y, MPFR_RNDN); 1509 /* As x > 0, the rounded value of x^y to nearest in precision p is equal 1510 to 0 iff x^y <= 2^(emin - 2). In particular, this does not depend on 1511 the precision p. Hence the following test. */ 1512 if (MPFR_IS_ZERO (z1) && MPFR_NOTZERO (z2)) 1513 { 1514 printf ("Error in bug20080820\n"); 1515 exit (1); 1516 } 1517 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); 1518 set_emin (emin); 1519 } 1520 1521 static void 1522 bug20110320 (void) 1523 { 1524 mpfr_exp_t emin; 1525 mpfr_t x, y, z1, z2; 1526 int inex; 1527 unsigned int flags; 1528 1529 emin = mpfr_get_emin (); 1530 mpfr_set_emin (11); 1531 mpfr_inits2 (2, x, y, z1, z2, (mpfr_ptr) 0); 1532 mpfr_set_ui_2exp (x, 1, 215, MPFR_RNDN); 1533 mpfr_set_ui (y, 1024, MPFR_RNDN); 1534 mpfr_clear_flags (); 1535 inex = mpfr_pow (z1, x, y, MPFR_RNDN); 1536 flags = __gmpfr_flags; 1537 mpfr_set_ui_2exp (z2, 1, 215*1024, MPFR_RNDN); 1538 if (inex != 0 || flags != 0 || ! mpfr_equal_p (z1, z2)) 1539 { 1540 printf ("Error in bug20110320\n"); 1541 printf ("Expected inex = 0, flags = 0, z = "); 1542 mpfr_dump (z2); 1543 printf ("Got inex = %d, flags = %u, z = ", inex, flags); 1544 mpfr_dump (z1); 1545 exit (1); 1546 } 1547 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); 1548 set_emin (emin); 1549 } 1550 1551 int 1552 main (int argc, char **argv) 1553 { 1554 mpfr_prec_t p; 1555 1556 tests_start_mpfr (); 1557 1558 bug20071127 (); 1559 special (); 1560 particular_cases (); 1561 check_pow_ui (); 1562 check_pow_si (); 1563 check_special_pow_si (); 1564 pow_si_long_min (); 1565 for (p = 2; p < 100; p++) 1566 check_inexact (p); 1567 underflows (); 1568 overflows (); 1569 overflows2 (); 1570 overflows3 (); 1571 x_near_one (); 1572 bug20071103 (); 1573 bug20071104 (); 1574 bug20071128 (); 1575 bug20071218 (); 1576 bug20080721 (); 1577 bug20080820 (); 1578 bug20110320 (); 1579 1580 test_generic (2, 100, 100); 1581 test_generic_ui (2, 100, 100); 1582 test_generic_si (2, 100, 100); 1583 1584 data_check ("data/pow275", mpfr_pow275, "mpfr_pow275"); 1585 1586 tests_end_mpfr (); 1587 return 0; 1588 } 1589