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