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