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