1 /* Test file for mpfr_root (also used for mpfr_rootn_ui). 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 #ifdef TF 24 # define TF_IS_MPFR_ROOT 0 25 #else 26 # define TF_IS_MPFR_ROOT 1 27 # define TF mpfr_root 28 # define _MPFR_NO_DEPRECATED_ROOT 29 #endif 30 31 #include "mpfr-test.h" 32 33 #include <time.h> 34 35 /* return the cpu time in seconds */ 36 static double 37 cputime (void) 38 { 39 return (double) clock () / (double) CLOCKS_PER_SEC; 40 } 41 42 #define DEFN(N) \ 43 static int root##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \ 44 { return TF (y, x, N, rnd); } \ 45 static int pow##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \ 46 { return mpfr_pow_ui (y, x, N, rnd); } 47 48 DEFN(2) 49 DEFN(3) 50 DEFN(4) 51 DEFN(5) 52 DEFN(17) 53 DEFN(120) 54 55 static void 56 special (void) 57 { 58 mpfr_t x, y; 59 int i; 60 61 mpfr_init (x); 62 mpfr_init (y); 63 64 /* root(NaN) = NaN */ 65 mpfr_set_nan (x); 66 i = TF (y, x, 17, MPFR_RNDN); 67 if (!mpfr_nan_p (y)) 68 { 69 printf ("Error: root(NaN,17) <> NaN\n"); 70 exit (1); 71 } 72 MPFR_ASSERTN (i == 0); 73 74 /* root(+Inf) = +Inf */ 75 mpfr_set_inf (x, 1); 76 i = TF (y, x, 42, MPFR_RNDN); 77 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) 78 { 79 printf ("Error: root(+Inf,42) <> +Inf\n"); 80 exit (1); 81 } 82 MPFR_ASSERTN (i == 0); 83 84 /* root(-Inf, 17) = -Inf */ 85 mpfr_set_inf (x, -1); 86 i = TF (y, x, 17, MPFR_RNDN); 87 if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0) 88 { 89 printf ("Error: root(-Inf,17) <> -Inf\n"); 90 exit (1); 91 } 92 MPFR_ASSERTN (i == 0); 93 94 /* root(-Inf, 42) = NaN */ 95 mpfr_set_inf (x, -1); 96 i = TF (y, x, 42, MPFR_RNDN); 97 if (!mpfr_nan_p (y)) 98 { 99 printf ("Error: root(-Inf,42) <> -Inf\n"); 100 exit (1); 101 } 102 MPFR_ASSERTN (i == 0); 103 104 /* root(+/-0, k) = +/-0, with the sign depending on TF. 105 * Before calling the function, we set y to NaN with the wrong sign, 106 * so that if the code of the function forgets to do something, this 107 * will be detected. 108 */ 109 mpfr_set_zero (x, 1); /* x is +0 */ 110 MPFR_SET_NAN (y); 111 MPFR_SET_NEG (y); 112 i = TF (y, x, 17, MPFR_RNDN); 113 if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y)) 114 { 115 printf ("Error: root(+0,17) <> +0\n"); 116 exit (1); 117 } 118 MPFR_ASSERTN (i == 0); 119 MPFR_SET_NAN (y); 120 MPFR_SET_NEG (y); 121 i = TF (y, x, 42, MPFR_RNDN); 122 if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y)) 123 { 124 printf ("Error: root(+0,42) <> +0\n"); 125 exit (1); 126 } 127 MPFR_ASSERTN (i == 0); 128 mpfr_set_zero (x, -1); /* x is -0 */ 129 MPFR_SET_NAN (y); 130 MPFR_SET_POS (y); 131 i = TF (y, x, 17, MPFR_RNDN); 132 if (MPFR_NOTZERO (y) || MPFR_IS_POS (y)) 133 { 134 printf ("Error: root(-0,17) <> -0\n"); 135 exit (1); 136 } 137 MPFR_ASSERTN (i == 0); 138 MPFR_SET_NAN (y); 139 if (TF_IS_MPFR_ROOT) 140 MPFR_SET_POS (y); 141 else 142 MPFR_SET_NEG (y); 143 i = TF (y, x, 42, MPFR_RNDN); 144 if (MPFR_NOTZERO (y) || 145 (TF_IS_MPFR_ROOT ? MPFR_IS_POS (y) : MPFR_IS_NEG (y))) 146 { 147 printf ("Error: root(-0,42) <> %c0\n", 148 TF_IS_MPFR_ROOT ? '-' : '+'); 149 exit (1); 150 } 151 MPFR_ASSERTN (i == 0); 152 153 mpfr_set_prec (x, 53); 154 mpfr_set_str (x, "8.39005285514734966412e-01", 10, MPFR_RNDN); 155 TF (x, x, 3, MPFR_RNDN); 156 if (mpfr_cmp_str1 (x, "9.43166207799662426048e-01")) 157 { 158 printf ("Error in root3 (1)\n"); 159 printf ("expected 9.43166207799662426048e-01\n"); 160 printf ("got "); 161 mpfr_dump (x); 162 exit (1); 163 } 164 165 mpfr_set_prec (x, 32); 166 mpfr_set_prec (y, 32); 167 mpfr_set_str_binary (x, "0.10000100001100101001001001011001"); 168 TF (x, x, 3, MPFR_RNDN); 169 mpfr_set_str_binary (y, "0.11001101011000100111000111111001"); 170 if (mpfr_cmp (x, y)) 171 { 172 printf ("Error in root3 (2)\n"); 173 exit (1); 174 } 175 176 mpfr_set_prec (x, 32); 177 mpfr_set_prec (y, 32); 178 mpfr_set_str_binary (x, "-0.1100001110110000010101011001011"); 179 TF (x, x, 3, MPFR_RNDD); 180 mpfr_set_str_binary (y, "-0.11101010000100100101000101011001"); 181 if (mpfr_cmp (x, y)) 182 { 183 printf ("Error in root3 (3)\n"); 184 exit (1); 185 } 186 187 mpfr_set_prec (x, 82); 188 mpfr_set_prec (y, 27); 189 mpfr_set_str_binary (x, "0.1010001111011101011011000111001011001101100011110110010011011011011010011001100101e-7"); 190 TF (y, x, 3, MPFR_RNDD); 191 mpfr_set_str_binary (x, "0.101011110001110001000100011E-2"); 192 if (mpfr_cmp (x, y)) 193 { 194 printf ("Error in root3 (4)\n"); 195 exit (1); 196 } 197 198 mpfr_set_prec (x, 204); 199 mpfr_set_prec (y, 38); 200 mpfr_set_str_binary (x, "0.101000000001101000000001100111111011111001110110100001111000100110100111001101100111110001110001011011010110010011100101111001111100001010010100111011101100000011011000101100010000000011000101001010001001E-5"); 201 TF (y, x, 3, MPFR_RNDD); 202 mpfr_set_str_binary (x, "0.10001001111010011011101000010110110010E-1"); 203 if (mpfr_cmp (x, y)) 204 { 205 printf ("Error in root3 (5)\n"); 206 exit (1); 207 } 208 209 /* Worst case found on 2006-11-25 */ 210 mpfr_set_prec (x, 53); 211 mpfr_set_prec (y, 53); 212 mpfr_set_str_binary (x, "1.0100001101101101001100110001001000000101001101100011E28"); 213 TF (y, x, 35, MPFR_RNDN); 214 mpfr_set_str_binary (x, "1.1100000010110101100011101011000010100001101100100011E0"); 215 if (mpfr_cmp (x, y)) 216 { 217 printf ("Error in rootn (y, x, 35, MPFR_RNDN) for\n" 218 "x = 1.0100001101101101001100110001001000000101001101100011E28\n" 219 "Expected "); 220 mpfr_dump (x); 221 printf ("Got "); 222 mpfr_dump (y); 223 exit (1); 224 } 225 /* Worst cases found on 2006-11-26 */ 226 mpfr_set_str_binary (x, "1.1111010011101110001111010110000101110000110110101100E17"); 227 TF (y, x, 36, MPFR_RNDD); 228 mpfr_set_str_binary (x, "1.0110100111010001101001010111001110010100111111000010E0"); 229 if (mpfr_cmp (x, y)) 230 { 231 printf ("Error in rootn (y, x, 36, MPFR_RNDD) for\n" 232 "x = 1.1111010011101110001111010110000101110000110110101100E17\n" 233 "Expected "); 234 mpfr_dump (x); 235 printf ("Got "); 236 mpfr_dump (y); 237 exit (1); 238 } 239 mpfr_set_str_binary (x, "1.1100011101101101100010110001000001110001111110010000E23"); 240 TF (y, x, 36, MPFR_RNDU); 241 mpfr_set_str_binary (x, "1.1001010100001110000110111111100011011101110011000100E0"); 242 if (mpfr_cmp (x, y)) 243 { 244 printf ("Error in rootn (y, x, 36, MPFR_RNDU) for\n" 245 "x = 1.1100011101101101100010110001000001110001111110010000E23\n" 246 "Expected "); 247 mpfr_dump (x); 248 printf ("Got "); 249 mpfr_dump (y); 250 exit (1); 251 } 252 253 /* Check for k = 1 */ 254 mpfr_set_ui (x, 17, MPFR_RNDN); 255 i = TF (y, x, 1, MPFR_RNDN); 256 if (mpfr_cmp_ui (x, 17) || i != 0) 257 { 258 printf ("Error in root for 17^(1/1)\n"); 259 exit (1); 260 } 261 262 mpfr_set_ui (x, 0, MPFR_RNDN); 263 i = TF (y, x, 0, MPFR_RNDN); 264 if (!MPFR_IS_NAN (y) || i != 0) 265 { 266 printf ("Error in root for (+0)^(1/0)\n"); 267 exit (1); 268 } 269 mpfr_neg (x, x, MPFR_RNDN); 270 i = TF (y, x, 0, MPFR_RNDN); 271 if (!MPFR_IS_NAN (y) || i != 0) 272 { 273 printf ("Error in root for (-0)^(1/0)\n"); 274 exit (1); 275 } 276 277 mpfr_set_ui (x, 1, MPFR_RNDN); 278 i = TF (y, x, 0, MPFR_RNDN); 279 if (!MPFR_IS_NAN (y) || i != 0) 280 { 281 printf ("Error in root for 1^(1/0)\n"); 282 exit (1); 283 } 284 285 /* Check for k==2 */ 286 mpfr_set_si (x, -17, MPFR_RNDD); 287 i = TF (y, x, 2, MPFR_RNDN); 288 if (!MPFR_IS_NAN (y) || i != 0) 289 { 290 printf ("Error in root for (-17)^(1/2)\n"); 291 exit (1); 292 } 293 294 mpfr_clear (x); 295 mpfr_clear (y); 296 } 297 298 /* https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=812779 299 * https://bugzilla.gnome.org/show_bug.cgi?id=756960 300 * is a GNOME Calculator bug (mpfr_root applied on a negative integer, 301 * which is converted to an unsigned integer), but the strange result 302 * is also due to a bug in MPFR. 303 */ 304 static void 305 bigint (void) 306 { 307 mpfr_t x, y; 308 309 mpfr_inits2 (64, x, y, (mpfr_ptr) 0); 310 311 mpfr_set_ui (x, 10, MPFR_RNDN); 312 if (sizeof (unsigned long) * CHAR_BIT == 64) 313 { 314 TF (x, x, ULONG_MAX, MPFR_RNDN); 315 mpfr_set_ui_2exp (y, 1, -63, MPFR_RNDN); 316 mpfr_add_ui (y, y, 1, MPFR_RNDN); 317 if (! mpfr_equal_p (x, y)) 318 { 319 printf ("Error in bigint for ULONG_MAX\n"); 320 printf ("Expected "); 321 mpfr_dump (y); 322 printf ("Got "); 323 mpfr_dump (x); 324 exit (1); 325 } 326 } 327 328 mpfr_set_ui (x, 10, MPFR_RNDN); 329 TF (x, x, 1234567890, MPFR_RNDN); 330 mpfr_set_str_binary (y, 331 "1.00000000000000000000000000001000000000101011000101000110010001"); 332 if (! mpfr_equal_p (x, y)) 333 { 334 printf ("Error in bigint for 1234567890\n"); 335 printf ("Expected "); 336 mpfr_dump (y); 337 printf ("Got "); 338 mpfr_dump (x); 339 exit (1); 340 } 341 342 mpfr_clears (x, y, (mpfr_ptr) 0); 343 } 344 345 #define TEST_FUNCTION TF 346 #define INTEGER_TYPE unsigned long 347 #define INT_RAND_FUNCTION() \ 348 (INTEGER_TYPE) (randlimb () & 1 ? randlimb () : randlimb () % 3 + 2) 349 #include "tgeneric_ui.c" 350 351 static void 352 exact_powers (unsigned long bmax, unsigned long kmax) 353 { 354 long b, k; 355 mpz_t z; 356 mpfr_t x, y; 357 int inex, neg; 358 359 mpz_init (z); 360 for (b = 2; b <= bmax; b++) 361 for (k = 1; k <= kmax; k++) 362 { 363 mpz_ui_pow_ui (z, b, k); 364 mpfr_init2 (x, mpz_sizeinbase (z, 2)); 365 mpfr_set_ui (x, b, MPFR_RNDN); 366 mpfr_pow_ui (x, x, k, MPFR_RNDN); 367 mpz_set_ui (z, b); 368 mpfr_init2 (y, mpz_sizeinbase (z, 2)); 369 for (neg = 0; neg <= 1; neg++) 370 { 371 inex = TF (y, x, k, MPFR_RNDN); 372 if (inex != 0) 373 { 374 printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k); 375 printf ("Expected inex=0, got %d\n", inex); 376 exit (1); 377 } 378 if (neg && (k & 1) == 0) 379 { 380 if (!MPFR_IS_NAN (y)) 381 { 382 printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k); 383 printf ("Expected y=NaN\n"); 384 printf ("Got "); 385 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN); 386 printf ("\n"); 387 exit (1); 388 } 389 } 390 else if (MPFR_IS_NAN (y) || mpfr_cmp_si (y, b) != 0) 391 { 392 printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k); 393 printf ("Expected y=%ld\n", b); 394 printf ("Got "); 395 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN); 396 printf ("\n"); 397 exit (1); 398 } 399 mpfr_neg (x, x, MPFR_RNDN); 400 b = -b; 401 } 402 mpfr_clear (x); 403 mpfr_clear (y); 404 } 405 mpz_clear (z); 406 } 407 408 /* Compare root(x,2^h) with pow(x,2^(-h)). */ 409 static void 410 cmp_pow (void) 411 { 412 mpfr_t x, y1, y2; 413 int h; 414 415 mpfr_inits2 (128, x, y1, y2, (mpfr_ptr) 0); 416 417 for (h = 1; h < sizeof (unsigned long) * CHAR_BIT; h++) 418 { 419 unsigned long k = (unsigned long) 1 << h; 420 int i; 421 422 for (i = 0; i < 10; i++) 423 { 424 mpfr_rnd_t rnd; 425 mpfr_flags_t flags1, flags2; 426 int inex1, inex2; 427 428 tests_default_random (x, 0, __gmpfr_emin, __gmpfr_emax, 1); 429 rnd = RND_RAND_NO_RNDF (); 430 mpfr_set_ui_2exp (y1, 1, -h, MPFR_RNDN); 431 mpfr_clear_flags (); 432 inex1 = mpfr_pow (y1, x, y1, rnd); 433 flags1 = __gmpfr_flags; 434 mpfr_clear_flags (); 435 inex2 = TF (y2, x, k, rnd); 436 flags2 = __gmpfr_flags; 437 if (!(mpfr_equal_p (y1, y2) && SAME_SIGN (inex1, inex2) && 438 flags1 == flags2)) 439 { 440 printf ("Error in cmp_pow on h=%d, i=%d, rnd=%s\n", 441 h, i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 442 printf ("x = "); 443 mpfr_dump (x); 444 printf ("pow = "); 445 mpfr_dump (y1); 446 printf ("with inex = %d, flags =", inex1); 447 flags_out (flags1); 448 printf ("root = "); 449 mpfr_dump (y2); 450 printf ("with inex = %d, flags =", inex2); 451 flags_out (flags2); 452 exit (1); 453 } 454 } 455 } 456 457 mpfr_clears (x, y1, y2, (mpfr_ptr) 0); 458 } 459 460 static void 461 bug20171214 (void) 462 { 463 mpfr_t x, y; 464 int inex; 465 466 mpfr_init2 (x, 805); 467 mpfr_init2 (y, 837); 468 mpfr_set_ui (x, 1, MPFR_RNDN); 469 inex = TF (y, x, 120, MPFR_RNDN); 470 MPFR_ASSERTN (inex == 0); 471 MPFR_ASSERTN (mpfr_cmp_ui (y, 1) == 0); 472 mpfr_set_si (x, -1, MPFR_RNDN); 473 inex = TF (y, x, 121, MPFR_RNDN); 474 MPFR_ASSERTN (inex == 0); 475 MPFR_ASSERTN (mpfr_cmp_si (y, -1) == 0); 476 mpfr_clear (x); 477 mpfr_clear (y); 478 } 479 480 int 481 main (int argc, char *argv[]) 482 { 483 mpfr_t x; 484 int r; 485 mpfr_prec_t p; 486 unsigned long k; 487 488 if (argc == 3) /* troot prec k */ 489 { 490 double st1, st2; 491 unsigned long k; 492 int l; 493 mpfr_t y; 494 p = strtoul (argv[1], NULL, 10); 495 k = strtoul (argv[2], NULL, 10); 496 mpfr_init2 (x, p); 497 mpfr_init2 (y, p); 498 mpfr_const_pi (y, MPFR_RNDN); 499 TF (x, y, k, MPFR_RNDN); /* to warm up cache */ 500 st1 = cputime (); 501 for (l = 0; cputime () - st1 < 1.0; l++) 502 TF (x, y, k, MPFR_RNDN); 503 st1 = (cputime () - st1) / l; 504 printf ("%-15s took %.2es\n", MAKE_STR(TF), st1); 505 506 /* compare with x^(1/k) = exp(1/k*log(x)) */ 507 /* first warm up cache */ 508 mpfr_swap (x, y); 509 mpfr_log (y, x, MPFR_RNDN); 510 mpfr_div_ui (y, y, k, MPFR_RNDN); 511 mpfr_exp (y, y, MPFR_RNDN); 512 513 st2 = cputime (); 514 for (l = 0; cputime () - st2 < 1.0; l++) 515 { 516 mpfr_log (y, x, MPFR_RNDN); 517 mpfr_div_ui (y, y, k, MPFR_RNDN); 518 mpfr_exp (y, y, MPFR_RNDN); 519 } 520 st2 = (cputime () - st2) / l; 521 printf ("exp(1/k*log(x)) took %.2es\n", st2); 522 523 mpfr_clear (x); 524 mpfr_clear (y); 525 return 0; 526 } 527 528 tests_start_mpfr (); 529 530 bug20171214 (); 531 exact_powers (3, 1000); 532 special (); 533 bigint (); 534 cmp_pow (); 535 536 mpfr_init (x); 537 538 for (p = MPFR_PREC_MIN; p < 100; p++) 539 { 540 mpfr_set_prec (x, p); 541 for (r = 0; r < MPFR_RND_MAX; r++) 542 { 543 mpfr_set_ui (x, 1, MPFR_RNDN); 544 k = 2 + randlimb () % 4; /* 2 <= k <= 5 */ 545 TF (x, x, k, (mpfr_rnd_t) r); 546 if (mpfr_cmp_ui (x, 1)) 547 { 548 printf ("Error in rootn[%lu] for x=1, rnd=%s\ngot ", 549 k, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 550 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 551 printf ("\n"); 552 exit (1); 553 } 554 mpfr_set_si (x, -1, MPFR_RNDN); 555 if (k % 2) 556 { 557 TF (x, x, k, (mpfr_rnd_t) r); 558 if (mpfr_cmp_si (x, -1)) 559 { 560 printf ("Error in rootn[%lu] for x=-1, rnd=%s\ngot ", 561 k, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 562 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 563 printf ("\n"); 564 exit (1); 565 } 566 } 567 568 if (p >= 5) 569 { 570 int i; 571 for (i = -12; i <= 12; i++) 572 { 573 mpfr_set_ui (x, 27, MPFR_RNDN); 574 mpfr_mul_2si (x, x, 3*i, MPFR_RNDN); 575 TF (x, x, 3, MPFR_RNDN); 576 if (mpfr_cmp_si_2exp (x, 3, i)) 577 { 578 printf ("Error in rootn[3] for " 579 "x = 27.0 * 2^(%d), rnd=%s\ngot ", 580 3*i, mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 581 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 582 printf ("\ninstead of 3 * 2^(%d)\n", i); 583 exit (1); 584 } 585 } 586 } 587 } 588 } 589 mpfr_clear (x); 590 591 test_generic_ui (MPFR_PREC_MIN, 200, 30); 592 593 bad_cases (root2, pow2, "rootn[2]", 0, -256, 255, 4, 128, 80, 40); 594 bad_cases (root3, pow3, "rootn[3]", 256, -256, 255, 4, 128, 200, 40); 595 bad_cases (root4, pow4, "rootn[4]", 0, -256, 255, 4, 128, 320, 40); 596 bad_cases (root5, pow5, "rootn[5]", 256, -256, 255, 4, 128, 440, 40); 597 bad_cases (root17, pow17, "rootn[17]", 256, -256, 255, 4, 128, 800, 40); 598 bad_cases (root120, pow120, "rootn[120]", 0, -256, 255, 4, 128, 800, 40); 599 600 tests_end_mpfr (); 601 return 0; 602 } 603