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 Free Software Foundation, Inc. 4 Contributed by the Arenaire and Cacao 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 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 197 static void 198 check_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 */ 226 static void 227 bug20071104 (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 259 static void 260 check_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) */ 299 static void 300 bug20080223 (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 321 static void 322 bug20080904 (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 351 int 352 main (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