1 /* Test file for mpfr_pow_z -- power function x^z with z a MPZ 2 3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 #include <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 31 static void 32 check_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 140 static void 141 check_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 (unsigned long) 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 = ", (unsigned long) p, n, 185 mpfr_print_rnd_mode (rnd)); 186 mpfr_dump (x); 187 printf ("Ypowsi(inex = %2d) = ", res1); mpfr_dump (y1); 188 printf ("Ypowz (inex = %2d) = ", res2); mpfr_dump (y2); 189 exit (1); 190 } 191 } 192 } /* for i */ 193 } /* for p */ 194 mpfr_clears (x, y1, y2, (mpfr_ptr) 0); 195 mpz_clear (z); 196 } 197 198 static void 199 check_regression (void) 200 { 201 mpfr_t x, y; 202 mpz_t z; 203 int res1, res2; 204 205 mpz_init_set_ui (z, 2026876995); 206 mpfr_init2 (x, 122); 207 mpfr_init2 (y, 122); 208 209 mpfr_set_str_binary (x, "0.10000010010000111101001110100101101010011110011100001111000001001101000110011001001001001011001011010110110110101000111011E1"); 210 res1 = mpfr_pow_z (y, x, z, MPFR_RNDU); 211 res2 = mpfr_pow_ui (x, x, 2026876995UL, MPFR_RNDU); 212 if (mpfr_cmp (x, y) || res1 != res2) 213 { 214 printf ("Regression (1) tested failed (%d=?%d)\n",res1, res2); 215 printf ("pow_ui: "); mpfr_dump (x); 216 printf ("pow_z: "); mpfr_dump (y); 217 218 exit (1); 219 } 220 221 mpfr_clear (x); 222 mpfr_clear (y); 223 mpz_clear (z); 224 } 225 226 /* Bug found by Kevin P. Rauch */ 227 static void 228 bug20071104 (void) 229 { 230 mpfr_t x, y; 231 mpz_t z; 232 int inex; 233 234 mpz_init_set_si (z, -2); 235 mpfr_inits2 (20, x, y, (mpfr_ptr) 0); 236 mpfr_set_ui (x, 0, MPFR_RNDN); 237 mpfr_nextbelow (x); /* x = -2^(emin-1) */ 238 mpfr_clear_flags (); 239 inex = mpfr_pow_z (y, x, z, MPFR_RNDN); 240 if (! mpfr_inf_p (y) || MPFR_SIGN (y) < 0) 241 { 242 printf ("Error in bug20071104: expected +Inf, got "); 243 mpfr_dump (y); 244 exit (1); 245 } 246 if (inex <= 0) 247 { 248 printf ("Error in bug20071104: bad ternary value (%d)\n", inex); 249 exit (1); 250 } 251 if (__gmpfr_flags != (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT)) 252 { 253 printf ("Error in bug20071104: bad flags (%u)\n", __gmpfr_flags); 254 exit (1); 255 } 256 mpfr_clears (x, y, (mpfr_ptr) 0); 257 mpz_clear (z); 258 } 259 260 static void 261 check_overflow (void) 262 { 263 mpfr_t a; 264 mpz_t z; 265 unsigned long n; 266 int res; 267 268 mpfr_init2 (a, 53); 269 270 mpfr_set_str_binary (a, "1E10"); 271 mpz_init_set_ui (z, ULONG_MAX); 272 res = mpfr_pow_z (a, a, z, MPFR_RNDN); 273 if (! MPFR_IS_INF (a) || MPFR_SIGN (a) < 0 || res <= 0) 274 { 275 printf ("Error for (1e10)^ULONG_MAX, expected +Inf,\ngot "); 276 mpfr_dump (a); 277 exit (1); 278 } 279 280 /* Bug in pow_z.c up to r5109: if x = y (same mpfr_t argument), the 281 input argument is negative and not a power of two, z is positive 282 and odd, an overflow or underflow occurs, and the temporary result 283 res is positive, then the result gets a wrong sign (positive 284 instead of negative). */ 285 mpfr_set_str_binary (a, "-1.1E10"); 286 n = (ULONG_MAX ^ (ULONG_MAX >> 1)) + 1; 287 mpz_set_ui (z, n); 288 res = mpfr_pow_z (a, a, z, MPFR_RNDN); 289 if (! MPFR_IS_INF (a) || MPFR_SIGN (a) > 0 || res >= 0) 290 { 291 printf ("Error for (-1e10)^%lu, expected -Inf,\ngot ", n); 292 mpfr_dump (a); 293 exit (1); 294 } 295 296 mpfr_clear (a); 297 mpz_clear (z); 298 } 299 300 /* bug reported by Carl Witty (32-bit architecture) */ 301 static void 302 bug20080223 (void) 303 { 304 mpfr_t a, exp, answer; 305 306 mpfr_init2 (a, 53); 307 mpfr_init2 (exp, 53); 308 mpfr_init2 (answer, 53); 309 310 mpfr_set_si (exp, -1073741824, MPFR_RNDN); 311 312 mpfr_set_str (a, "1.999999999", 10, MPFR_RNDN); 313 /* a = 562949953139837/2^48 */ 314 mpfr_pow (answer, a, exp, MPFR_RNDN); 315 mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823"); 316 MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0); 317 318 mpfr_clear (a); 319 mpfr_clear (exp); 320 mpfr_clear (answer); 321 } 322 323 static void 324 bug20080904 (void) 325 { 326 mpz_t exp; 327 mpfr_t a, answer; 328 mpfr_exp_t emin_default; 329 330 mpz_init (exp); 331 mpfr_init2 (a, 70); 332 mpfr_init2 (answer, 70); 333 334 emin_default = mpfr_get_emin (); 335 mpfr_set_emin (MPFR_EMIN_MIN); 336 337 mpz_set_str (exp, "-4eb92f8c7b7bf81e", 16); 338 mpfr_set_str_binary (a, "1.110000101110100110100011111000011110111101000011111001111001010011100"); 339 340 mpfr_pow_z (answer, a, exp, MPFR_RNDN); 341 /* The correct result is near 2^(-2^62), so it underflows when 342 MPFR_EMIN_MIN > -2^62 (i.e. with 32 and 64 bits machines). */ 343 mpfr_set_str (a, "AA500C0D7A69275DBp-4632850503556296886", 16, MPFR_RNDN); 344 if (! mpfr_equal_p (answer, a)) 345 { 346 printf ("Error in bug20080904:\n"); 347 printf ("Expected "); 348 mpfr_out_str (stdout, 16, 0, a, MPFR_RNDN); 349 putchar ('\n'); 350 printf ("Got "); 351 mpfr_out_str (stdout, 16, 0, answer, MPFR_RNDN); 352 putchar ('\n'); 353 exit (1); 354 } 355 356 mpfr_set_emin (emin_default); 357 358 mpz_clear (exp); 359 mpfr_clear (a); 360 mpfr_clear (answer); 361 } 362 363 int 364 main (void) 365 { 366 tests_start_mpfr (); 367 368 check_special (); 369 370 check_integer (2, 163, 100); 371 check_regression (); 372 bug20071104 (); 373 bug20080223 (); 374 bug20080904 (); 375 check_overflow (); 376 377 tests_end_mpfr (); 378 return 0; 379 } 380