1/* Test file for mpfr_pow_z -- power function x^z with z a MPZ 2 3Copyright 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 <stdlib.h> 24#include <float.h> 25#include <math.h> 26 27#include "mpfr-test.h" 28 29#define ERROR(str) do { printf("Error for "str"\n"); exit (1); } while (0) 30 31static void 32check_special (void) 33{ 34 mpfr_t x, y; 35 mpz_t z; 36 int res; 37 38 mpfr_init (x); 39 mpfr_init (y); 40 mpz_init (z); 41 42 /* x^0 = 1 except for NAN */ 43 mpfr_set_ui (x, 23, MPFR_RNDN); 44 mpz_set_ui (z, 0); 45 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 46 if (res != 0 || mpfr_cmp_ui (y, 1) != 0) 47 ERROR ("23^0"); 48 mpfr_set_nan (x); 49 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 50 if (res != 0 || mpfr_nan_p (y) || mpfr_cmp_si (y, 1) != 0) 51 ERROR ("NAN^0"); 52 mpfr_set_inf (x, 1); 53 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 54 if (res != 0 || mpfr_cmp_ui (y, 1) != 0) 55 ERROR ("INF^0"); 56 57 /* sINF^N = INF if s==1 or n even if N > 0*/ 58 mpz_set_ui (z, 42); 59 mpfr_set_inf (x, 1); 60 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 61 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0) 62 ERROR ("INF^42"); 63 mpfr_set_inf (x, -1); 64 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 65 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0) 66 ERROR ("-INF^42"); 67 mpz_set_ui (z, 17); 68 mpfr_set_inf (x, 1); 69 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 70 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) <= 0) 71 ERROR ("INF^17"); 72 mpfr_set_inf (x, -1); 73 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 74 if (res != 0 || mpfr_inf_p (y) == 0 || mpfr_sgn (y) >= 0) 75 ERROR ("-INF^17"); 76 77 mpz_set_si (z, -42); 78 mpfr_set_inf (x, 1); 79 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 80 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) 81 ERROR ("INF^-42"); 82 mpfr_set_inf (x, -1); 83 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 84 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) 85 ERROR ("-INF^-42"); 86 mpz_set_si (z, -17); 87 mpfr_set_inf (x, 1); 88 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 89 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) 90 ERROR ("INF^-17"); 91 mpfr_set_inf (x, -1); 92 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 93 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) >= 0) 94 ERROR ("-INF^-17"); 95 96 /* s0^N = +0 if s==+ or n even if N > 0*/ 97 mpz_set_ui (z, 42); 98 MPFR_SET_ZERO (x); MPFR_SET_POS (x); 99 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 100 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) 101 ERROR ("+0^42"); 102 MPFR_SET_NEG (x); 103 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 104 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) 105 ERROR ("-0^42"); 106 mpz_set_ui (z, 17); 107 MPFR_SET_POS (x); 108 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 109 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) <= 0) 110 ERROR ("+0^17"); 111 MPFR_SET_NEG (x); 112 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 113 if (res != 0 || mpfr_zero_p (y) == 0 || MPFR_SIGN (y) >= 0) 114 ERROR ("-0^17"); 115 116 mpz_set_si (z, -42); 117 MPFR_SET_ZERO (x); MPFR_SET_POS (x); 118 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 119 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) <= 0) 120 ERROR ("+0^-42"); 121 MPFR_SET_NEG (x); 122 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 123 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) <= 0) 124 ERROR ("-0^-42"); 125 mpz_set_si (z, -17); 126 MPFR_SET_POS (x); 127 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 128 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) <= 0) 129 ERROR ("+0^-17"); 130 MPFR_SET_NEG (x); 131 res = mpfr_pow_z (y, x, z, MPFR_RNDN); 132 if (res != 0 || mpfr_inf_p (y) == 0 || MPFR_SIGN (y) >= 0) 133 ERROR ("-0^-17"); 134 135 mpz_clear (z); 136 mpfr_clear (y); 137 mpfr_clear (x); 138} 139 140static void 141check_integer (mpfr_prec_t begin, mpfr_prec_t end, unsigned long max) 142{ 143 mpfr_t x, y1, y2; 144 mpz_t z; 145 unsigned long i, n; 146 mpfr_prec_t p; 147 int res1, res2; 148 mpfr_rnd_t rnd; 149 150 mpfr_inits2 (begin, x, y1, y2, (mpfr_ptr) 0); 151 mpz_init (z); 152 for (p = begin ; p < end ; p+=4) 153 { 154 mpfr_set_prec (x, p); 155 mpfr_set_prec (y1, p); 156 mpfr_set_prec (y2, p); 157 for (i = 0 ; i < max ; i++) 158 { 159 mpz_urandomb (z, RANDS, GMP_NUMB_BITS); 160 if ((i & 1) != 0) 161 mpz_neg (z, z); 162 mpfr_urandomb (x, RANDS); 163 mpfr_mul_2ui (x, x, 1, MPFR_RNDN); /* 0 <= x < 2 */ 164 rnd = RND_RAND (); 165 if (mpz_fits_slong_p (z)) 166 { 167 n = mpz_get_si (z); 168 /* printf ("New test for x=%ld\nCheck Pow_si\n", n); */ 169 res1 = mpfr_pow_si (y1, x, n, rnd); 170 /* printf ("Check pow_z\n"); */ 171 res2 = mpfr_pow_z (y2, x, z, rnd); 172 if (mpfr_cmp (y1, y2) != 0) 173 { 174 printf ("Error for p = %lu, z = %lu, rnd = %s and x = ", 175 p, n, mpfr_print_rnd_mode (rnd)); 176 mpfr_dump (x); 177 printf ("Ypowsi = "); mpfr_dump (y1); 178 printf ("Ypowz = "); mpfr_dump (y2); 179 exit (1); 180 } 181 if (res1 != res2) 182 { 183 printf ("Wrong inexact flags for p = %lu, z = %lu, rnd = %s" 184 " and x = ", p, n, mpfr_print_rnd_mode (rnd)); 185 mpfr_dump (x); 186 printf ("Ypowsi(inex = %2d) = ", res1); mpfr_dump (y1); 187 printf ("Ypowz (inex = %2d) = ", res2); mpfr_dump (y2); 188 exit (1); 189 } 190 } 191 } /* for i */ 192 } /* for p */ 193 mpfr_clears (x, y1, y2, (mpfr_ptr) 0); 194 mpz_clear (z); 195} 196 197static void 198check_regression (void) 199{ 200 mpfr_t x, y; 201 mpz_t z; 202 int res1, res2; 203 204 mpz_init_set_ui (z, 2026876995); 205 mpfr_init2 (x, 122); 206 mpfr_init2 (y, 122); 207 208 mpfr_set_str_binary (x, "0.10000010010000111101001110100101101010011110011100001111000001001101000110011001001001001011001011010110110110101000111011E1"); 209 res1 = mpfr_pow_z (y, x, z, MPFR_RNDU); 210 res2 = mpfr_pow_ui (x, x, 2026876995UL, MPFR_RNDU); 211 if (mpfr_cmp (x, y) || res1 != res2) 212 { 213 printf ("Regression (1) tested failed (%d=?%d)\n",res1, res2); 214 printf ("pow_ui: "); mpfr_dump (x); 215 printf ("pow_z: "); mpfr_dump (y); 216 217 exit (1); 218 } 219 220 mpfr_clear (x); 221 mpfr_clear (y); 222 mpz_clear (z); 223} 224 225/* Bug found by Kevin P. Rauch */ 226static void 227bug20071104 (void) 228{ 229 mpfr_t x, y; 230 mpz_t z; 231 int inex; 232 233 mpz_init_set_si (z, -2); 234 mpfr_inits2 (20, x, y, (mpfr_ptr) 0); 235 mpfr_set_ui (x, 0, MPFR_RNDN); 236 mpfr_nextbelow (x); /* x = -2^(emin-1) */ 237 mpfr_clear_flags (); 238 inex = mpfr_pow_z (y, x, z, MPFR_RNDN); 239 if (! mpfr_inf_p (y) || MPFR_SIGN (y) < 0) 240 { 241 printf ("Error in bug20071104: expected +Inf, got "); 242 mpfr_dump (y); 243 exit (1); 244 } 245 if (inex <= 0) 246 { 247 printf ("Error in bug20071104: bad ternary value (%d)\n", inex); 248 exit (1); 249 } 250 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 251 { 252 printf ("Error in bug20071104: bad flags (%u)\n", __gmpfr_flags); 253 exit (1); 254 } 255 mpfr_clears (x, y, (mpfr_ptr) 0); 256 mpz_clear (z); 257} 258 259static void 260check_overflow (void) 261{ 262 mpfr_t a; 263 mpz_t z; 264 unsigned long n; 265 int res; 266 267 mpfr_init2 (a, 53); 268 269 mpfr_set_str_binary (a, "1E10"); 270 mpz_init_set_ui (z, ULONG_MAX); 271 res = mpfr_pow_z (a, a, z, MPFR_RNDN); 272 if (!MPFR_IS_INF (a) || MPFR_SIGN (a) < 0) 273 { 274 printf ("Error for (1e10)^ULONG_MAX\n"); 275 exit (1); 276 } 277 278 /* Bug in pow_z.c up to r5109: if x = y (same mpfr_t argument), the 279 input argument is negative and not a power of two, z is positive 280 and odd, an overflow or underflow occurs, and the temporary result 281 res is positive, then the result gets a wrong sign (positive 282 instead of negative). */ 283 mpfr_set_str_binary (a, "-1.1E10"); 284 n = (ULONG_MAX ^ (ULONG_MAX >> 1)) + 1; 285 mpz_set_ui (z, n); 286 res = mpfr_pow_z (a, a, z, MPFR_RNDN); 287 if (!MPFR_IS_INF (a) || MPFR_SIGN (a) > 0) 288 { 289 printf ("Error for (-1e10)^%lu, expected -Inf,\ngot ", n); 290 mpfr_dump (a); 291 exit (1); 292 } 293 294 mpfr_clear (a); 295 mpz_clear (z); 296} 297 298/* bug reported by Carl Witty (32-bit architecture) */ 299static void 300bug20080223 (void) 301{ 302 mpfr_t a, exp, answer; 303 304 mpfr_init2 (a, 53); 305 mpfr_init2 (exp, 53); 306 mpfr_init2 (answer, 53); 307 308 mpfr_set_si (exp, -1073741824, MPFR_RNDN); 309 310 mpfr_set_str (a, "1.999999999", 10, MPFR_RNDN); 311 /* a = 562949953139837/2^48 */ 312 mpfr_pow (answer, a, exp, MPFR_RNDN); 313 mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823"); 314 MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0); 315 316 mpfr_clear (a); 317 mpfr_clear (exp); 318 mpfr_clear (answer); 319} 320 321static void 322bug20080904 (void) 323{ 324 mpz_t exp; 325 mpfr_t a, answer; 326 mpfr_exp_t emin_default; 327 328 mpz_init (exp); 329 mpfr_init2 (a, 70); 330 mpfr_init2 (answer, 70); 331 332 emin_default = mpfr_get_emin (); 333 mpfr_set_emin (MPFR_EMIN_MIN); 334 335 mpz_set_str (exp, "-4eb92f8c7b7bf81e", 16); 336 mpfr_set_str_binary (a, "1.110000101110100110100011111000011110111101000011111001111001010011100"); 337 338 mpfr_pow_z (answer, a, exp, MPFR_RNDN); 339 /* The correct result is near 2^(-2^62), so it underflows when 340 MPFR_EMIN_MIN > -2^62 (i.e. with 32 and 64 bits machines). */ 341 mpfr_set_str (a, "AA500C0D7A69275DBp-4632850503556296886", 16, MPFR_RNDN); 342 MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0); 343 344 mpfr_set_emin (emin_default); 345 346 mpz_clear (exp); 347 mpfr_clear (a); 348 mpfr_clear (answer); 349} 350 351int 352main (void) 353{ 354 tests_start_mpfr (); 355 356 check_special (); 357 358 check_integer (2, 163, 100); 359 check_regression (); 360 bug20071104 (); 361 bug20080223 (); 362 bug20080904 (); 363 check_overflow (); 364 365 tests_end_mpfr (); 366 return 0; 367} 368