1/* Test file for mpfr_pow, mpfr_pow_ui and mpfr_pow_si. 2 3Copyright 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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#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 32static int 33test_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 74static void 75check_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 260static void 261check_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, 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, 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, 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, - (LONG_MIN + 1))); 344 } 345 346 mpfr_clear (x); 347} 348 349static void 350check_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 384static void 385pow_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 413static void 414check_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 477static void 478special (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 705static void 706particular_cases (void) 707{ 708 mpfr_t t[11], r; 709 static const char *name[11] = { 710 "NaN", "+inf", "-inf", "+0", "-0", "+1", "-1", "+2", "-2", "+0.5", "-0.5"}; 711 int i, j; 712 int error = 0; 713 714 for (i = 0; i < 11; i++) 715 mpfr_init2 (t[i], 2); 716 mpfr_init2 (r, 6); 717 718 mpfr_set_nan (t[0]); 719 mpfr_set_inf (t[1], 1); 720 mpfr_set_ui (t[3], 0, MPFR_RNDN); 721 mpfr_set_ui (t[5], 1, MPFR_RNDN); 722 mpfr_set_ui (t[7], 2, MPFR_RNDN); 723 mpfr_div_2ui (t[9], t[5], 1, MPFR_RNDN); 724 for (i = 1; i < 11; i += 2) 725 mpfr_neg (t[i+1], t[i], MPFR_RNDN); 726 727 for (i = 0; i < 11; i++) 728 for (j = 0; j < 11; j++) 729 { 730 double d; 731 int p; 732 static int q[11][11] = { 733 /* NaN +inf -inf +0 -0 +1 -1 +2 -2 +0.5 -0.5 */ 734 /* NaN */ { 0, 0, 0, 128, 128, 0, 0, 0, 0, 0, 0 }, 735 /* +inf */ { 0, 1, 2, 128, 128, 1, 2, 1, 2, 1, 2 }, 736 /* -inf */ { 0, 1, 2, 128, 128, -1, -2, 1, 2, 1, 2 }, 737 /* +0 */ { 0, 2, 1, 128, 128, 2, 1, 2, 1, 2, 1 }, 738 /* -0 */ { 0, 2, 1, 128, 128, -2, -1, 2, 1, 2, 1 }, 739 /* +1 */ {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128 }, 740 /* -1 */ { 0, 128, 128, 128, 128,-128,-128, 128, 128, 0, 0 }, 741 /* +2 */ { 0, 1, 2, 128, 128, 256, 64, 512, 32, 180, 90 }, 742 /* -2 */ { 0, 1, 2, 128, 128,-256, -64, 512, 32, 0, 0 }, 743 /* +0.5 */ { 0, 2, 1, 128, 128, 64, 256, 32, 512, 90, 180 }, 744 /* -0.5 */ { 0, 2, 1, 128, 128, -64,-256, 32, 512, 0, 0 } 745 }; 746 test_pow (r, t[i], t[j], MPFR_RNDN); 747 p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 : 748 mpfr_cmp_ui (r, 0) == 0 ? 2 : 749 (d = mpfr_get_d (r, MPFR_RNDN), (int) (ABS(d) * 128.0)); 750 if (p != 0 && MPFR_SIGN(r) < 0) 751 p = -p; 752 if (p != q[i][j]) 753 { 754 printf ("Error in mpfr_pow for particular case (%s)^(%s) (%d,%d):\n" 755 "got %d instead of %d\n", name[i], name[j], i,j,p, q[i][j]); 756 mpfr_dump (r); 757 error = 1; 758 } 759 } 760 761 for (i = 0; i < 11; i++) 762 mpfr_clear (t[i]); 763 mpfr_clear (r); 764 765 if (error) 766 exit (1); 767} 768 769static void 770underflows (void) 771{ 772 mpfr_t x, y, z; 773 int err = 0; 774 int inexact; 775 int i; 776 mpfr_exp_t emin; 777 778 mpfr_init2 (x, 64); 779 mpfr_init2 (y, 64); 780 781 mpfr_set_ui (x, 1, MPFR_RNDN); 782 mpfr_set_exp (x, mpfr_get_emin()); 783 784 for (i = 3; i < 10; i++) 785 { 786 mpfr_set_ui (y, i, MPFR_RNDN); 787 mpfr_div_2ui (y, y, 1, MPFR_RNDN); 788 test_pow (y, x, y, MPFR_RNDN); 789 if (!MPFR_IS_FP(y) || mpfr_cmp_ui (y, 0)) 790 { 791 printf ("Error in mpfr_pow for "); 792 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 793 printf (" ^ (%d/2)\nGot ", i); 794 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 795 printf (" instead of 0.\n"); 796 exit (1); 797 } 798 } 799 800 mpfr_init2 (z, 55); 801 mpfr_set_str (x, "0.110011010011101001110001110100010000110111101E0", 802 2, MPFR_RNDN); 803 mpfr_set_str (y, "0.101110010011111001011010100011011100111110011E40", 804 2, MPFR_RNDN); 805 mpfr_clear_flags (); 806 inexact = mpfr_pow (z, x, y, MPFR_RNDU); 807 if (!mpfr_underflow_p ()) 808 { 809 printf ("Underflow flag is not set for special underflow test.\n"); 810 err = 1; 811 } 812 if (inexact <= 0) 813 { 814 printf ("Ternary value is wrong for special underflow test.\n"); 815 err = 1; 816 } 817 mpfr_set_ui (x, 0, MPFR_RNDN); 818 mpfr_nextabove (x); 819 if (mpfr_cmp (x, z) != 0) 820 { 821 printf ("Wrong value for special underflow test.\nGot "); 822 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN); 823 printf ("\ninstead of "); 824 mpfr_out_str (stdout, 2, 2, x, MPFR_RNDN); 825 printf ("\n"); 826 err = 1; 827 } 828 if (err) 829 exit (1); 830 831 /* MPFR currently (2006-08-19) segfaults on the following code (and 832 possibly makes other programs crash due to the lack of memory), 833 because y is converted into an mpz_t, and the required precision 834 is too high. */ 835 mpfr_set_prec (x, 2); 836 mpfr_set_prec (y, 2); 837 mpfr_set_prec (z, 12); 838 mpfr_set_ui_2exp (x, 3, -2, MPFR_RNDN); 839 mpfr_set_ui_2exp (y, 1, mpfr_get_emax () - 1, MPFR_RNDN); 840 mpfr_clear_flags (); 841 mpfr_pow (z, x, y, MPFR_RNDN); 842 if (!mpfr_underflow_p () || MPFR_NOTZERO (z)) 843 { 844 printf ("Underflow test with large y fails.\n"); 845 exit (1); 846 } 847 848 emin = mpfr_get_emin (); 849 mpfr_set_emin (-256); 850 mpfr_set_prec (x, 2); 851 mpfr_set_prec (y, 2); 852 mpfr_set_prec (z, 12); 853 mpfr_set_ui_2exp (x, 3, -2, MPFR_RNDN); 854 mpfr_set_ui_2exp (y, 1, 38, MPFR_RNDN); 855 mpfr_clear_flags (); 856 inexact = mpfr_pow (z, x, y, MPFR_RNDN); 857 if (!mpfr_underflow_p () || MPFR_NOTZERO (z) || inexact >= 0) 858 { 859 printf ("Bad underflow detection for 0.75^(2^38). Obtained:\n" 860 "Underflow flag... %-3s (should be 'yes')\n" 861 "Zero result...... %-3s (should be 'yes')\n" 862 "Inexact value.... %-3d (should be negative)\n", 863 mpfr_underflow_p () ? "yes" : "no", 864 MPFR_IS_ZERO (z) ? "yes" : "no", inexact); 865 exit (1); 866 } 867 mpfr_set_emin (emin); 868 869 emin = mpfr_get_emin (); 870 mpfr_set_emin (-256); 871 mpfr_set_prec (x, 2); 872 mpfr_set_prec (y, 40); 873 mpfr_set_prec (z, 12); 874 mpfr_set_ui_2exp (x, 3, -1, MPFR_RNDN); 875 mpfr_set_si_2exp (y, -1, 38, MPFR_RNDN); 876 for (i = 0; i < 4; i++) 877 { 878 if (i == 2) 879 mpfr_neg (x, x, MPFR_RNDN); 880 mpfr_clear_flags (); 881 inexact = mpfr_pow (z, x, y, MPFR_RNDN); 882 if (!mpfr_underflow_p () || MPFR_NOTZERO (z) || 883 (i == 3 ? (inexact <= 0) : (inexact >= 0))) 884 { 885 printf ("Bad underflow detection for ("); 886 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); 887 printf (")^(-2^38-%d). Obtained:\n" 888 "Overflow flag.... %-3s (should be 'no')\n" 889 "Underflow flag... %-3s (should be 'yes')\n" 890 "Zero result...... %-3s (should be 'yes')\n" 891 "Inexact value.... %-3d (should be %s)\n", i, 892 mpfr_overflow_p () ? "yes" : "no", 893 mpfr_underflow_p () ? "yes" : "no", 894 MPFR_IS_ZERO (z) ? "yes" : "no", inexact, 895 i == 3 ? "positive" : "negative"); 896 exit (1); 897 } 898 inexact = mpfr_sub_ui (y, y, 1, MPFR_RNDN); 899 MPFR_ASSERTN (inexact == 0); 900 } 901 mpfr_set_emin (emin); 902 903 mpfr_clears (x, y, z, (mpfr_ptr) 0); 904} 905 906static void 907overflows (void) 908{ 909 mpfr_t a, b; 910 911 /* bug found by Ming J. Tsai <mingjt@delvron.us>, 4 Oct 2003 */ 912 913 mpfr_init_set_str (a, "5.1e32", 10, MPFR_RNDN); 914 mpfr_init (b); 915 916 test_pow (b, a, a, MPFR_RNDN); 917 if (!(mpfr_inf_p (b) && mpfr_sgn (b) > 0)) 918 { 919 printf ("Error for a^a for a=5.1e32\n"); 920 printf ("Expected +Inf, got "); 921 mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN); 922 printf ("\n"); 923 exit (1); 924 } 925 926 mpfr_clear(a); 927 mpfr_clear(b); 928} 929 930static void 931overflows2 (void) 932{ 933 mpfr_t x, y, z; 934 mpfr_exp_t emin, emax; 935 int e; 936 937 /* x^y in reduced exponent range, where x = 2^b and y is not an integer 938 (so that mpfr_pow_z is not used). */ 939 940 emin = mpfr_get_emin (); 941 emax = mpfr_get_emax (); 942 set_emin (-128); 943 944 mpfr_inits2 (16, x, y, z, (mpfr_ptr) 0); 945 946 mpfr_set_si_2exp (x, 1, -64, MPFR_RNDN); /* 2^(-64) */ 947 mpfr_set_si_2exp (y, -1, -1, MPFR_RNDN); /* -0.5 */ 948 for (e = 2; e <= 32; e += 17) 949 { 950 set_emax (e); 951 mpfr_clear_flags (); 952 mpfr_pow (z, x, y, MPFR_RNDN); 953 if (MPFR_IS_NEG (z) || ! mpfr_inf_p (z)) 954 { 955 printf ("Error in overflows2 (e = %d): expected +Inf, got ", e); 956 mpfr_dump (z); 957 exit (1); 958 } 959 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 960 { 961 printf ("Error in overflows2 (e = %d): bad flags (%u)\n", 962 e, __gmpfr_flags); 963 exit (1); 964 } 965 } 966 967 mpfr_clears (x, y, z, (mpfr_ptr) 0); 968 969 set_emin (emin); 970 set_emax (emax); 971} 972 973static void 974overflows3 (void) 975{ 976 /* x^y where x = 2^b, y is not an integer (so that mpfr_pow_z is not used) 977 and b * y = emax in the extended exponent range. If emax is divisible 978 by 3, we choose x = 2^(-2*emax/3) and y = -3/2. 979 Test also with nextbelow(x). */ 980 981 if (MPFR_EMAX_MAX % 3 == 0) 982 { 983 mpfr_t x, y, z, t; 984 mpfr_exp_t emin, emax; 985 unsigned int flags; 986 int i; 987 988 emin = mpfr_get_emin (); 989 emax = mpfr_get_emax (); 990 set_emin (MPFR_EMIN_MIN); 991 set_emax (MPFR_EMAX_MAX); 992 993 mpfr_inits2 (16, x, y, z, t, (mpfr_ptr) 0); 994 995 mpfr_set_si_2exp (x, 1, -2 * (MPFR_EMAX_MAX / 3), MPFR_RNDN); 996 for (i = 0; i <= 1; i++) 997 { 998 mpfr_set_si_2exp (y, -3, -1, MPFR_RNDN); 999 mpfr_clear_flags (); 1000 mpfr_pow (z, x, y, MPFR_RNDN); 1001 if (MPFR_IS_NEG (z) || ! mpfr_inf_p (z)) 1002 { 1003 printf ("Error in overflows3 (RNDN, i = %d): expected +Inf," 1004 " got ", i); 1005 mpfr_dump (z); 1006 exit (1); 1007 } 1008 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 1009 { 1010 printf ("Error in overflows3 (RNDN, i = %d): bad flags (%u)\n", 1011 i, __gmpfr_flags); 1012 exit (1); 1013 } 1014 1015 mpfr_clear_flags (); 1016 mpfr_pow (z, x, y, MPFR_RNDZ); 1017 flags = __gmpfr_flags; 1018 mpfr_set (t, z, MPFR_RNDN); 1019 mpfr_nextabove (t); 1020 if (MPFR_IS_NEG (z) || mpfr_inf_p (z) || ! mpfr_inf_p (t)) 1021 { 1022 printf ("Error in overflows3 (RNDZ, i = %d):\nexpected ", i); 1023 mpfr_set_inf (t, 1); 1024 mpfr_nextbelow (t); 1025 mpfr_dump (t); 1026 printf ("got "); 1027 mpfr_dump (z); 1028 exit (1); 1029 } 1030 if (flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 1031 { 1032 printf ("Error in overflows3 (RNDZ, i = %d): bad flags (%u)\n", 1033 i, flags); 1034 exit (1); 1035 } 1036 mpfr_nextbelow (x); 1037 } 1038 1039 mpfr_clears (x, y, z, t, (mpfr_ptr) 0); 1040 1041 set_emin (emin); 1042 set_emax (emax); 1043 } 1044} 1045 1046static void 1047x_near_one (void) 1048{ 1049 mpfr_t x, y, z; 1050 int inex; 1051 1052 mpfr_init2 (x, 32); 1053 mpfr_init2 (y, 4); 1054 mpfr_init2 (z, 33); 1055 1056 mpfr_set_ui (x, 1, MPFR_RNDN); 1057 mpfr_nextbelow (x); 1058 mpfr_set_ui_2exp (y, 11, -2, MPFR_RNDN); 1059 inex = mpfr_pow (z, x, y, MPFR_RNDN); 1060 if (mpfr_cmp_str (z, "0.111111111111111111111111111111011E0", 2, MPFR_RNDN) 1061 || inex <= 0) 1062 { 1063 printf ("Failure in x_near_one, got inex = %d and\nz = ", inex); 1064 mpfr_dump (z); 1065 } 1066 1067 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1068} 1069 1070static int 1071mpfr_pow275 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) 1072{ 1073 mpfr_t z; 1074 int inex; 1075 1076 mpfr_init2 (z, 4); 1077 mpfr_set_ui_2exp (z, 11, -2, MPFR_RNDN); 1078 inex = mpfr_pow (y, x, z, MPFR_RNDN); 1079 mpfr_clear (z); 1080 return inex; 1081} 1082 1083/* Bug found by Kevin P. Rauch */ 1084static void 1085bug20071103 (void) 1086{ 1087 mpfr_t x, y, z; 1088 mpfr_exp_t emin, emax; 1089 1090 emin = mpfr_get_emin (); 1091 emax = mpfr_get_emax (); 1092 mpfr_set_emin (-1000000); 1093 mpfr_set_emax ( 1000000); 1094 1095 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0); 1096 mpfr_set_si_2exp (x, -3, -1, MPFR_RNDN); /* x = -1.5 */ 1097 mpfr_set_str (y, "-0.ffffffffffffffff", 16, MPFR_RNDN); 1098 mpfr_set_exp (y, mpfr_get_emax ()); 1099 mpfr_clear_flags (); 1100 mpfr_pow (z, x, y, MPFR_RNDN); 1101 MPFR_ASSERTN (mpfr_zero_p (z) && MPFR_SIGN (z) > 0 && 1102 __gmpfr_flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT)); 1103 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1104 1105 set_emin (emin); 1106 set_emax (emax); 1107} 1108 1109/* Bug found by Kevin P. Rauch */ 1110static void 1111bug20071104 (void) 1112{ 1113 mpfr_t x, y, z; 1114 mpfr_exp_t emin, emax; 1115 int inex; 1116 1117 emin = mpfr_get_emin (); 1118 emax = mpfr_get_emax (); 1119 mpfr_set_emin (-1000000); 1120 mpfr_set_emax ( 1000000); 1121 1122 mpfr_inits2 (20, x, y, z, (mpfr_ptr) 0); 1123 mpfr_set_ui (x, 0, MPFR_RNDN); 1124 mpfr_nextbelow (x); /* x = -2^(emin-1) */ 1125 mpfr_set_si (y, -2, MPFR_RNDN); /* y = -2 */ 1126 mpfr_clear_flags (); 1127 inex = mpfr_pow (z, x, y, MPFR_RNDN); 1128 if (! mpfr_inf_p (z) || MPFR_SIGN (z) < 0) 1129 { 1130 printf ("Error in bug20071104: expected +Inf, got "); 1131 mpfr_dump (z); 1132 exit (1); 1133 } 1134 if (inex <= 0) 1135 { 1136 printf ("Error in bug20071104: bad ternary value (%d)\n", inex); 1137 exit (1); 1138 } 1139 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 1140 { 1141 printf ("Error in bug20071104: bad flags (%u)\n", __gmpfr_flags); 1142 exit (1); 1143 } 1144 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1145 1146 set_emin (emin); 1147 set_emax (emax); 1148} 1149 1150/* Bug found by Kevin P. Rauch */ 1151static void 1152bug20071127 (void) 1153{ 1154 mpfr_t x, y, z; 1155 int i, tern; 1156 mpfr_exp_t emin, emax; 1157 1158 emin = mpfr_get_emin (); 1159 emax = mpfr_get_emax (); 1160 mpfr_set_emin (-1000000); 1161 mpfr_set_emax ( 1000000); 1162 1163 mpfr_init2 (x, 128); 1164 mpfr_init2 (y, 128); 1165 mpfr_init2 (z, 128); 1166 1167 mpfr_set_str (x, "0.80000000000000000000000000000001", 16, MPFR_RNDN); 1168 1169 for (i = 1; i < 9; i *= 2) 1170 { 1171 mpfr_set_str (y, "8000000000000000", 16, MPFR_RNDN); 1172 mpfr_add_si (y, y, i, MPFR_RNDN); 1173 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1174 MPFR_ASSERTN(mpfr_zero_p (z) && MPFR_IS_POS(z)); 1175 } 1176 1177 mpfr_clear (x); 1178 mpfr_clear (y); 1179 mpfr_clear (z); 1180 1181 mpfr_set_emin (emin); 1182 mpfr_set_emax (emax); 1183} 1184 1185/* Bug found by Kevin P. Rauch */ 1186static void 1187bug20071128 (void) 1188{ 1189 mpfr_t max_val, x, y, z; 1190 int i, tern; 1191 mpfr_exp_t emin, emax; 1192 1193 emin = mpfr_get_emin (); 1194 emax = mpfr_get_emax (); 1195 mpfr_set_emin (-1000000); 1196 mpfr_set_emax ( 1000000); 1197 1198 mpfr_init2 (max_val, 64); 1199 mpfr_init2 (x, 64); 1200 mpfr_init2 (y, 64); 1201 mpfr_init2 (z, 64); 1202 1203 mpfr_set_str (max_val, "0.ffffffffffffffff", 16, MPFR_RNDN); 1204 mpfr_set_exp (max_val, mpfr_get_emax ()); 1205 1206 mpfr_neg (x, max_val, MPFR_RNDN); 1207 1208 /* on 64-bit machines */ 1209 for (i = 41; i < 45; i++) 1210 { 1211 mpfr_set_si_2exp (y, -1, i, MPFR_RNDN); 1212 mpfr_add_si (y, y, 1, MPFR_RNDN); 1213 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1214 MPFR_ASSERTN(mpfr_zero_p (z) && MPFR_SIGN(z) < 0); 1215 } 1216 1217 /* on 32-bit machines */ 1218 for (i = 9; i < 13; i++) 1219 { 1220 mpfr_set_si_2exp (y, -1, i, MPFR_RNDN); 1221 mpfr_add_si (y, y, 1, MPFR_RNDN); 1222 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1223 MPFR_ASSERTN(mpfr_zero_p (z) && MPFR_SIGN(z) < 0); 1224 } 1225 1226 mpfr_clear (x); 1227 mpfr_clear (y); 1228 mpfr_clear (z); 1229 mpfr_clear (max_val); 1230 1231 mpfr_set_emin (emin); 1232 mpfr_set_emax (emax); 1233} 1234 1235/* Bug found by Kevin P. Rauch */ 1236static void 1237bug20071218 (void) 1238{ 1239 mpfr_t x, y, z, t; 1240 int tern; 1241 1242 mpfr_inits2 (64, x, y, z, t, (mpfr_ptr) 0); 1243 mpfr_set_str (x, "0x.80000000000002P-1023", 0, MPFR_RNDN); 1244 mpfr_set_str (y, "100000.000000002", 16, MPFR_RNDN); 1245 mpfr_set_ui (t, 0, MPFR_RNDN); 1246 mpfr_nextabove (t); 1247 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1248 if (mpfr_cmp0 (z, t) != 0) 1249 { 1250 printf ("Error in bug20071218 (1): Expected\n"); 1251 mpfr_dump (t); 1252 printf ("Got\n"); 1253 mpfr_dump (z); 1254 exit (1); 1255 } 1256 if (tern <= 0) 1257 { 1258 printf ("Error in bug20071218 (1): bad ternary value" 1259 " (%d instead of positive)\n", tern); 1260 exit (1); 1261 } 1262 mpfr_mul_2ui (y, y, 32, MPFR_RNDN); 1263 tern = mpfr_pow (z, x, y, MPFR_RNDN); 1264 if (MPFR_NOTZERO (z) || MPFR_IS_NEG (z)) 1265 { 1266 printf ("Error in bug20071218 (2): expected 0, got\n"); 1267 mpfr_dump (z); 1268 exit (1); 1269 } 1270 if (tern >= 0) 1271 { 1272 printf ("Error in bug20071218 (2): bad ternary value" 1273 " (%d instead of negative)\n", tern); 1274 exit (1); 1275 } 1276 mpfr_clears (x, y, z, t, (mpfr_ptr) 0); 1277} 1278 1279/* With revision 5429, this gives: 1280 * pow.c:43: assertion failed: !mpfr_integer_p (y) 1281 * This is fixed in revision 5432. 1282 */ 1283static void 1284bug20080721 (void) 1285{ 1286 mpfr_t x, y, z, t[2]; 1287 int inex; 1288 int rnd; 1289 int err = 0; 1290 1291 /* Note: input values have been chosen in a way to select the 1292 * general case. If mpfr_pow is modified, in particular line 1293 * if (y_is_integer && (MPFR_GET_EXP (y) <= 256)) 1294 * make sure that this test still does what we want. 1295 */ 1296 mpfr_inits2 (4913, x, y, (mpfr_ptr) 0); 1297 mpfr_inits2 (8, z, t[0], t[1], (mpfr_ptr) 0); 1298 mpfr_set_si (x, -1, MPFR_RNDN); 1299 mpfr_nextbelow (x); 1300 mpfr_set_ui_2exp (y, 1, mpfr_get_prec (y) - 1, MPFR_RNDN); 1301 inex = mpfr_add_ui (y, y, 1, MPFR_RNDN); 1302 MPFR_ASSERTN (inex == 0); 1303 mpfr_set_str_binary (t[0], "-0.10101101e2"); 1304 mpfr_set_str_binary (t[1], "-0.10101110e2"); 1305 RND_LOOP (rnd) 1306 { 1307 int i, inex0; 1308 1309 i = (rnd == MPFR_RNDN || rnd == MPFR_RNDD || rnd == MPFR_RNDA); 1310 inex0 = i ? -1 : 1; 1311 mpfr_clear_flags (); 1312 inex = mpfr_pow (z, x, y, (mpfr_rnd_t) rnd); 1313 if (__gmpfr_flags != MPFR_FLAGS_INEXACT || ! SAME_SIGN (inex, inex0) 1314 || MPFR_IS_NAN (z) || mpfr_cmp (z, t[i]) != 0) 1315 { 1316 unsigned int flags = __gmpfr_flags; 1317 1318 printf ("Error in bug20080721 with %s\n", 1319 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 1320 printf ("expected "); 1321 mpfr_out_str (stdout, 2, 0, t[i], MPFR_RNDN); 1322 printf (", inex = %d, flags = %u\n", inex0, 1323 (unsigned int) MPFR_FLAGS_INEXACT); 1324 printf ("got "); 1325 mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN); 1326 printf (", inex = %d, flags = %u\n", inex, flags); 1327 err = 1; 1328 } 1329 } 1330 mpfr_clears (x, y, z, t[0], t[1], (mpfr_ptr) 0); 1331 if (err) 1332 exit (1); 1333} 1334 1335/* The following test fails in r5552 (32-bit and 64-bit). This is due to: 1336 * mpfr_log (t, absx, MPFR_RNDU); 1337 * mpfr_mul (t, y, t, MPFR_RNDU); 1338 * in pow.c, that is supposed to compute an upper bound on exp(y*ln|x|), 1339 * but this is incorrect if y is negative. 1340 */ 1341static void 1342bug20080820 (void) 1343{ 1344 mpfr_exp_t emin; 1345 mpfr_t x, y, z1, z2; 1346 1347 emin = mpfr_get_emin (); 1348 mpfr_set_emin (MPFR_EMIN_MIN); 1349 mpfr_init2 (x, 80); 1350 mpfr_init2 (y, sizeof (mpfr_exp_t) * CHAR_BIT + 32); 1351 mpfr_init2 (z1, 2); 1352 mpfr_init2 (z2, 80); 1353 mpfr_set_ui (x, 2, MPFR_RNDN); 1354 mpfr_nextbelow (x); 1355 mpfr_set_exp_t (y, mpfr_get_emin () - 2, MPFR_RNDN); 1356 mpfr_nextabove (y); 1357 mpfr_pow (z1, x, y, MPFR_RNDN); 1358 mpfr_pow (z2, x, y, MPFR_RNDN); 1359 /* As x > 0, the rounded value of x^y to nearest in precision p is equal 1360 to 0 iff x^y <= 2^(emin - 2). In particular, this does not depend on 1361 the precision p. Hence the following test. */ 1362 if (MPFR_IS_ZERO (z1) && MPFR_NOTZERO (z2)) 1363 { 1364 printf ("Error in bug20080820\n"); 1365 exit (1); 1366 } 1367 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); 1368 set_emin (emin); 1369} 1370 1371static void 1372bug20110320 (void) 1373{ 1374 mpfr_exp_t emin; 1375 mpfr_t x, y, z1, z2; 1376 int inex; 1377 unsigned int flags; 1378 1379 emin = mpfr_get_emin (); 1380 mpfr_set_emin (11); 1381 mpfr_inits2 (2, x, y, z1, z2, (mpfr_ptr) 0); 1382 mpfr_set_ui_2exp (x, 1, 215, MPFR_RNDN); 1383 mpfr_set_ui (y, 1024, MPFR_RNDN); 1384 mpfr_clear_flags (); 1385 inex = mpfr_pow (z1, x, y, MPFR_RNDN); 1386 flags = __gmpfr_flags; 1387 mpfr_set_ui_2exp (z2, 1, 215*1024, MPFR_RNDN); 1388 if (inex != 0 || flags != 0 || ! mpfr_equal_p (z1, z2)) 1389 { 1390 printf ("Error in bug20110320\n"); 1391 printf ("Expected inex = 0, flags = 0, z = "); 1392 mpfr_dump (z2); 1393 printf ("Got inex = %d, flags = %u, z = ", inex, flags); 1394 mpfr_dump (z1); 1395 exit (1); 1396 } 1397 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); 1398 set_emin (emin); 1399} 1400 1401int 1402main (int argc, char **argv) 1403{ 1404 mpfr_prec_t p; 1405 1406 tests_start_mpfr (); 1407 1408 bug20071127 (); 1409 special (); 1410 particular_cases (); 1411 check_pow_ui (); 1412 check_pow_si (); 1413 check_special_pow_si (); 1414 pow_si_long_min (); 1415 for (p = 2; p < 100; p++) 1416 check_inexact (p); 1417 underflows (); 1418 overflows (); 1419 overflows2 (); 1420 overflows3 (); 1421 x_near_one (); 1422 bug20071103 (); 1423 bug20071104 (); 1424 bug20071128 (); 1425 bug20071218 (); 1426 bug20080721 (); 1427 bug20080820 (); 1428 bug20110320 (); 1429 1430 test_generic (2, 100, 100); 1431 test_generic_ui (2, 100, 100); 1432 test_generic_si (2, 100, 100); 1433 1434 data_check ("data/pow275", mpfr_pow275, "mpfr_pow275"); 1435 1436 tests_end_mpfr (); 1437 return 0; 1438} 1439