1 /* Test file for mpfr_compound_si. 2 3 Copyright 2021-2023 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 TEST_FUNCTION mpfr_compound_si 26 #define INTEGER_TYPE long 27 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS) 28 #define test_generic_ui test_generic_si 29 #include "tgeneric_ui.c" 30 31 /* Special cases from IEEE 754-2019 */ 32 static void 33 check_ieee754 (void) 34 { 35 mpfr_t x, y; 36 long i; 37 long t[] = { 0, 1, 2, 3, 17, LONG_MAX-1, LONG_MAX }; 38 int j; 39 mpfr_prec_t prec = 2; /* we need at least 2 so that 3/4 is exact */ 40 41 mpfr_init2 (x, prec); 42 mpfr_init2 (y, prec); 43 44 /* compound(x,n) = NaN for x < -1, and set invalid exception */ 45 for (i = 0; i < numberof(t); i++) 46 for (j = 0; j < 2; j++) 47 { 48 const char *s; 49 50 mpfr_clear_nanflag (); 51 if (j == 0) 52 { 53 mpfr_set_si (x, -2, MPFR_RNDN); 54 s = "-2"; 55 } 56 else 57 { 58 mpfr_set_inf (x, -1); 59 s = "-Inf"; 60 } 61 mpfr_compound_si (y, x, t[i], MPFR_RNDN); 62 if (!mpfr_nan_p (y)) 63 { 64 printf ("Error, compound(%s,%ld) should give NaN\n", s, t[i]); 65 printf ("got "); mpfr_dump (y); 66 exit (1); 67 } 68 if (!mpfr_nanflag_p ()) 69 { 70 printf ("Error, compound(%s,%ld) should raise invalid flag\n", 71 s, t[i]); 72 exit (1); 73 } 74 } 75 76 /* compound(x,0) = 1 for x >= -1 or x = NaN */ 77 for (i = -2; i <= 2; i++) 78 { 79 if (i == -2) 80 mpfr_set_nan (x); 81 else if (i == 2) 82 mpfr_set_inf (x, 1); 83 else 84 mpfr_set_si (x, i, MPFR_RNDN); 85 mpfr_compound_si (y, x, 0, MPFR_RNDN); 86 if (mpfr_cmp_ui (y, 1) != 0) 87 { 88 printf ("Error, compound(x,0) should give 1 on\nx = "); 89 mpfr_dump (x); 90 printf ("got "); mpfr_dump (y); 91 exit (1); 92 } 93 } 94 95 /* compound(-1,n) = +Inf for n < 0, and raise divide-by-zero flag */ 96 mpfr_clear_divby0 (); 97 mpfr_set_si (x, -1, MPFR_RNDN); 98 mpfr_compound_si (y, x, -1, MPFR_RNDN); 99 if (!mpfr_inf_p (y) || MPFR_SIGN(y) < 0) 100 { 101 printf ("Error, compound(-1,-1) should give +Inf\n"); 102 printf ("got "); mpfr_dump (y); 103 exit (1); 104 } 105 if (!mpfr_divby0_p ()) 106 { 107 printf ("Error, compound(-1,-1) should raise divide-by-zero flag\n"); 108 exit (1); 109 } 110 111 /* compound(-1,n) = +0 for n > 0 */ 112 mpfr_set_si (x, -1, MPFR_RNDN); 113 mpfr_compound_si (y, x, 1, MPFR_RNDN); 114 if (!mpfr_zero_p (y) || MPFR_SIGN(y) < 0) 115 { 116 printf ("Error, compound(-1,1) should give +0\n"); 117 printf ("got "); mpfr_dump (y); 118 exit (1); 119 } 120 121 /* compound(+/-0,n) = 1 */ 122 for (i = -1; i <= 1; i++) 123 { 124 mpfr_set_zero (x, -1); 125 mpfr_compound_si (y, x, i, MPFR_RNDN); 126 if (mpfr_cmp_ui (y, 1) != 0) 127 { 128 printf ("Error1, compound(x,%ld) should give 1\non x = ", i); 129 mpfr_dump (x); 130 printf ("got "); mpfr_dump (y); 131 exit (1); 132 } 133 mpfr_set_zero (x, +1); 134 mpfr_compound_si (y, x, i, MPFR_RNDN); 135 if (mpfr_cmp_ui (y, 1) != 0) 136 { 137 printf ("Error, compound(x,%ld) should give 1\non x = ", i); 138 mpfr_dump (x); 139 printf ("got "); mpfr_dump (y); 140 exit (1); 141 } 142 } 143 144 /* compound(+Inf,n) = +Inf for n > 0 */ 145 mpfr_set_inf (x, 1); 146 mpfr_compound_si (y, x, 1, MPFR_RNDN); 147 if (!mpfr_inf_p (y) || MPFR_SIGN(y) < 0) 148 { 149 printf ("Error, compound(+Inf,1) should give +Inf\n"); 150 printf ("got "); mpfr_dump (y); 151 exit (1); 152 } 153 154 /* compound(+Inf,n) = +0 for n < 0 */ 155 mpfr_set_inf (x, 1); 156 mpfr_compound_si (y, x, -1, MPFR_RNDN); 157 if (!mpfr_zero_p (y) || MPFR_SIGN(y) < 0) 158 { 159 printf ("Error, compound(+Inf,-1) should give +0\n"); 160 printf ("got "); mpfr_dump (y); 161 exit (1); 162 } 163 164 /* compound(NaN,n) = NaN for n <> 0 */ 165 mpfr_set_nan (x); 166 mpfr_compound_si (y, x, -1, MPFR_RNDN); 167 if (!mpfr_nan_p (y)) 168 { 169 printf ("Error, compound(NaN,-1) should give NaN\n"); 170 printf ("got "); mpfr_dump (y); 171 exit (1); 172 } 173 mpfr_compound_si (y, x, +1, MPFR_RNDN); 174 if (!mpfr_nan_p (y)) 175 { 176 printf ("Error, compound(NaN,+1) should give NaN\n"); 177 printf ("got "); mpfr_dump (y); 178 exit (1); 179 } 180 181 /* hard-coded test: x is the 32-bit nearest approximation of 17/42 */ 182 mpfr_set_prec (x, 32); 183 mpfr_set_prec (y, 32); 184 mpfr_set_ui_2exp (x, 3476878287UL, -33, MPFR_RNDN); 185 mpfr_compound_si (y, x, 12, MPFR_RNDN); 186 mpfr_set_ui_2exp (x, 1981447393UL, -25, MPFR_RNDN); 187 if (!mpfr_equal_p (y, x)) 188 { 189 printf ("Error for compound(3476878287/2^33,12)\n"); 190 printf ("expected "); mpfr_dump (x); 191 printf ("got "); mpfr_dump (y); 192 exit (1); 193 } 194 195 /* test for negative n */ 196 i = -1; 197 while (1) 198 { 199 /* i has the form -(2^k-1) */ 200 mpfr_set_si_2exp (x, -1, -1, MPFR_RNDN); /* x = -0.5 */ 201 mpfr_compound_si (y, x, i, MPFR_RNDN); 202 mpfr_set_ui_2exp (x, 1, -i, MPFR_RNDN); 203 if (!mpfr_equal_p (y, x)) 204 { 205 printf ("Error for compound(-0.5,%ld)\n", i); 206 printf ("expected "); mpfr_dump (x); 207 printf ("got "); mpfr_dump (y); 208 exit (1); 209 } 210 if (i == -2147483647) /* largest possible value on 32-bit machine */ 211 break; 212 i = 2 * i - 1; 213 } 214 215 /* The "#if" makes sure that 64-bit constants are supported, avoiding 216 a compilation failure. The "if" makes sure that the constant is 217 representable in a long (this would not be the case with 32-bit 218 unsigned long and 64-bit limb). */ 219 #if GMP_NUMB_BITS >= 64 || MPFR_PREC_BITS >= 64 220 if (4994322635099777669 <= LONG_MAX) 221 { 222 i = -4994322635099777669; 223 mpfr_set_ui (x, 1, MPFR_RNDN); 224 mpfr_compound_si (y, x, i, MPFR_RNDN); 225 mpfr_set_si (x, 1, MPFR_RNDN); 226 mpfr_mul_2si (x, x, i, MPFR_RNDN); 227 if (!mpfr_equal_p (y, x)) 228 { 229 printf ("Error for compound(1,%ld)\n", i); 230 printf ("expected "); mpfr_dump (x); 231 printf ("got "); mpfr_dump (y); 232 exit (1); 233 } 234 } 235 #endif 236 237 mpfr_clear (x); 238 mpfr_clear (y); 239 } 240 241 /* Failure with mpfr_compound_si from 2021-02-15 due to 242 incorrect underflow detection. */ 243 static void 244 bug_20230206 (void) 245 { 246 if (MPFR_PREC_MIN == 1) 247 { 248 mpfr_t x, y1, y2; 249 int inex1, inex2; 250 mpfr_flags_t flags1, flags2; 251 #if MPFR_PREC_BITS >= 64 252 mpfr_exp_t emin; 253 #endif 254 255 mpfr_inits2 (1, x, y1, y2, (mpfr_ptr) 0); 256 mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN); /* x = 1/2 */ 257 258 /* This first test is useful mainly for a 32-bit mpfr_exp_t type 259 (no failure with a 64-bit mpfr_exp_t type since the underflow 260 threshold in the extended exponent range is much lower). */ 261 262 mpfr_set_ui_2exp (y1, 1, -1072124363, MPFR_RNDN); 263 inex1 = -1; 264 flags1 = MPFR_FLAGS_INEXACT; 265 mpfr_clear_flags (); 266 /* -1832808704 ~= -2^30 / log2(3/2) */ 267 inex2 = mpfr_compound_si (y2, x, -1832808704, MPFR_RNDN); 268 flags2 = __gmpfr_flags; 269 if (!(mpfr_equal_p (y1, y2) && 270 SAME_SIGN (inex1, inex2) && 271 flags1 == flags2)) 272 { 273 printf ("Error in bug_20230206 (1):\n"); 274 printf ("Expected "); 275 mpfr_dump (y1); 276 printf (" with inex = %d, flags =", inex1); 277 flags_out (flags1); 278 printf ("Got "); 279 mpfr_dump (y2); 280 printf (" with inex = %d, flags =", inex2); 281 flags_out (flags2); 282 exit (1); 283 } 284 285 /* This second test is for a 64-bit mpfr_exp_t type 286 (it is disabled with a 32-bit mpfr_exp_t type). */ 287 288 /* The "#if" makes sure that 64-bit constants are supported, avoiding 289 a compilation failure. The "if" makes sure that the constant is 290 representable in a long (this would not be the case with 32-bit 291 unsigned long and 64-bit limb). It also ensures that mpfr_exp_t 292 has at least 64 bits. */ 293 #if MPFR_PREC_BITS >= 64 294 emin = mpfr_get_emin (); 295 set_emin (MPFR_EMIN_MIN); 296 mpfr_set_ui_2exp (y1, 1, -4611686018427366846, MPFR_RNDN); 297 inex1 = 1; 298 flags1 = MPFR_FLAGS_INEXACT; 299 mpfr_clear_flags (); 300 /* -7883729320669216768 ~= -2^62 / log2(3/2) */ 301 inex2 = mpfr_compound_si (y2, x, -7883729320669216768, MPFR_RNDN); 302 flags2 = __gmpfr_flags; 303 if (!(mpfr_equal_p (y1, y2) && 304 SAME_SIGN (inex1, inex2) && 305 flags1 == flags2)) 306 { 307 printf ("Error in bug_20230206 (2):\n"); 308 printf ("Expected "); 309 mpfr_dump (y1); 310 printf (" with inex = %d, flags =", inex1); 311 flags_out (flags1); 312 printf ("Got "); 313 mpfr_dump (y2); 314 printf (" with inex = %d, flags =", inex2); 315 flags_out (flags2); 316 exit (1); 317 } 318 set_emin (emin); 319 #endif 320 321 mpfr_clears (x, y1, y2, (mpfr_ptr) 0); 322 } 323 } 324 325 /* Reported by Patrick Pelissier on 2023-02-11 for the master branch 326 (tgeneric_ui.c with GMP_CHECK_RANDOMIZE=1412991715). 327 On a 32-bit host, one gets Inf (overflow) instead of 0.1E1071805703. 328 */ 329 static void 330 bug_20230211 (void) 331 { 332 mpfr_t x, y1, y2; 333 int inex1, inex2; 334 mpfr_flags_t flags1, flags2; 335 336 mpfr_inits2 (1, x, y1, y2, (mpfr_ptr) 0); 337 mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN); /* x = 1/2 */ 338 mpfr_set_ui_2exp (y1, 1, 1071805702, MPFR_RNDN); 339 inex1 = 1; 340 flags1 = MPFR_FLAGS_INEXACT; 341 mpfr_clear_flags (); 342 inex2 = mpfr_compound_si (y2, x, 1832263949, MPFR_RNDN); 343 flags2 = __gmpfr_flags; 344 if (!(mpfr_equal_p (y1, y2) && 345 SAME_SIGN (inex1, inex2) && 346 flags1 == flags2)) 347 { 348 printf ("Error in bug_20230211:\n"); 349 printf ("Expected "); 350 mpfr_dump (y1); 351 printf (" with inex = %d, flags =", inex1); 352 flags_out (flags1); 353 printf ("Got "); 354 mpfr_dump (y2); 355 printf (" with inex = %d, flags =", inex2); 356 flags_out (flags2); 357 exit (1); 358 } 359 mpfr_clears (x, y1, y2, (mpfr_ptr) 0); 360 } 361 362 /* Integer overflow with compound.c d04caeae04c6a83276916c4fbac1fe9b0cec3c8b 363 (2023-02-23) or 952fb0f5cc2df1fffde3eb54c462fdae5f123ea6 in the 4.2 branch 364 on "n * (kx - 1) + 1". Note: if the only effect is just a random value, 365 this probably doesn't affect the result (one might enter the "if" while 366 one shouldn't, but the real check is done inside the "if"). This test 367 fails if -fsanitize=undefined -fno-sanitize-recover is used or if the 368 processor emits a signal in case of integer overflow. 369 This test has been made obsolete by the "kx < ex" condition 370 in 2cb3123891dd46fe0258d4aec7f8655b8ec69aaf (master branch) 371 or f5cb40571bc3d1559f05b230cf4ffecaf0952852 (4.2 branch). */ 372 static void 373 bug_20230517 (void) 374 { 375 mpfr_exp_t old_emax; 376 mpfr_t x; 377 378 old_emax = mpfr_get_emax (); 379 set_emax (MPFR_EMAX_MAX); 380 381 mpfr_init2 (x, 123456); 382 mpfr_set_ui (x, 65536, MPFR_RNDN); 383 mpfr_nextabove (x); 384 mpfr_compound_si (x, x, LONG_MAX >> 16, MPFR_RNDN); 385 mpfr_clear (x); 386 387 set_emax (old_emax); 388 } 389 390 /* Inverse function on non-special cases... 391 One has x = (1+y)^n with y > -1 and x > 0. Thus y = x^(1/n) - 1. 392 The inverse function is useful 393 - to build and check hard-to-round cases (see bad_cases() in tests.c); 394 - to test the behavior close to the overflow and underflow thresholds. 395 The case x = 0 actually needs to be handled as it may occur with 396 bad_cases() due to rounding. 397 */ 398 static int 399 inv_compound (mpfr_ptr y, mpfr_srcptr x, long n, mpfr_rnd_t rnd_mode) 400 { 401 mpfr_t t; 402 int inexact; 403 mpfr_prec_t precy, prect; 404 MPFR_ZIV_DECL (loop); 405 MPFR_SAVE_EXPO_DECL (expo); 406 407 MPFR_ASSERTN (n != 0); 408 409 if (MPFR_UNLIKELY (MPFR_IS_ZERO (x))) 410 { 411 if (n > 0) 412 return mpfr_set_si (y, -1, rnd_mode); 413 else 414 { 415 MPFR_SET_INF (y); 416 MPFR_SET_POS (y); 417 MPFR_RET (0); 418 } 419 } 420 421 MPFR_SAVE_EXPO_MARK (expo); 422 423 if (mpfr_equal_p (x, __gmpfr_one)) 424 { 425 MPFR_SAVE_EXPO_FREE (expo); 426 mpfr_set_zero (y, 1); 427 MPFR_RET (0); 428 } 429 430 precy = MPFR_GET_PREC (y); 431 prect = precy + 20; 432 mpfr_init2 (t, prect); 433 434 MPFR_ZIV_INIT (loop, prect); 435 for (;;) 436 { 437 mpfr_exp_t expt1, expt2, err; 438 unsigned int inext; 439 440 if (mpfr_rootn_si (t, x, n, MPFR_RNDN) == 0) 441 { 442 /* With a huge t, this case would yield inext != 0 and a 443 MPFR_CAN_ROUND failure until a huge precision is reached 444 (as the result is very close to an exact point). Fortunately, 445 since t is exact, we can obtain the correctly rounded result 446 by doing the second operation to the target precision directly. 447 */ 448 inexact = mpfr_sub_ui (y, t, 1, rnd_mode); 449 goto end; 450 } 451 expt1 = MPFR_GET_EXP (t); 452 /* |error| <= 2^(expt1-prect-1) */ 453 inext = mpfr_sub_ui (t, t, 1, MPFR_RNDN); 454 if (MPFR_UNLIKELY (MPFR_IS_ZERO (t))) 455 goto cont; /* cannot round yet */ 456 expt2 = MPFR_GET_EXP (t); 457 err = 1; 458 if (expt2 < expt1) 459 err += expt1 - expt2; 460 /* |error(rootn)| <= 2^(err+expt2-prect-2) 461 and if mpfr_sub_ui is inexact: 462 |error| <= 2^(err+expt2-prect-2) + 2^(expt2-prect-1) 463 <= (2^(err-1) + 1) * 2^(expt2-prect-1) 464 <= 2^((err+1)+expt2-prect-2) */ 465 if (inext) 466 err++; 467 /* |error| <= 2^(err+expt2-prect-2) */ 468 if (MPFR_CAN_ROUND (t, prect + 2 - err, precy, rnd_mode)) 469 break; 470 471 cont: 472 MPFR_ZIV_NEXT (loop, prect); 473 mpfr_set_prec (t, prect); 474 } 475 476 inexact = mpfr_set (y, t, rnd_mode); 477 478 end: 479 MPFR_ZIV_FREE (loop); 480 mpfr_clear (t); 481 MPFR_SAVE_EXPO_FREE (expo); 482 return mpfr_check_range (y, inexact, rnd_mode); 483 } 484 485 #define DEFN(N) \ 486 static int mpfr_compound##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) \ 487 { return mpfr_compound_si (y, x, N, r); } \ 488 static int inv_compound##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) \ 489 { return inv_compound (y, x, N, r); } 490 491 DEFN(2) 492 DEFN(3) 493 DEFN(4) 494 DEFN(5) 495 DEFN(17) 496 DEFN(120) 497 498 #define TEST_FUNCTION mpfr_compound2 499 #define test_generic test_generic_compound2 500 #include "tgeneric.c" 501 502 #define TEST_FUNCTION mpfr_compound3 503 #define test_generic test_generic_compound3 504 #include "tgeneric.c" 505 506 #define TEST_FUNCTION mpfr_compound4 507 #define test_generic test_generic_compound4 508 #include "tgeneric.c" 509 510 #define TEST_FUNCTION mpfr_compound5 511 #define test_generic test_generic_compound5 512 #include "tgeneric.c" 513 514 #define TEST_FUNCTION mpfr_compound17 515 #define test_generic test_generic_compound17 516 #include "tgeneric.c" 517 518 #define TEST_FUNCTION mpfr_compound120 519 #define test_generic test_generic_compound120 520 #include "tgeneric.c" 521 522 int 523 main (void) 524 { 525 tests_start_mpfr (); 526 527 check_ieee754 (); 528 bug_20230206 (); 529 bug_20230211 (); 530 bug_20230517 (); 531 532 test_generic_si (MPFR_PREC_MIN, 100, 100); 533 534 test_generic_compound2 (MPFR_PREC_MIN, 100, 100); 535 test_generic_compound3 (MPFR_PREC_MIN, 100, 100); 536 test_generic_compound4 (MPFR_PREC_MIN, 100, 100); 537 test_generic_compound5 (MPFR_PREC_MIN, 100, 100); 538 test_generic_compound17 (MPFR_PREC_MIN, 100, 100); 539 test_generic_compound120 (MPFR_PREC_MIN, 100, 100); 540 541 /* Note: For small n, we need a psup high enough to avoid too many 542 "f exact while f^(-1) inexact" occurrences in bad_cases(). */ 543 bad_cases (mpfr_compound2, inv_compound2, "mpfr_compound2", 544 0, -256, 255, 4, 128, 240, 40); 545 bad_cases (mpfr_compound3, inv_compound3, "mpfr_compound3", 546 0, -256, 255, 4, 128, 120, 40); 547 bad_cases (mpfr_compound4, inv_compound4, "mpfr_compound4", 548 0, -256, 255, 4, 128, 80, 40); 549 bad_cases (mpfr_compound5, inv_compound5, "mpfr_compound5", 550 0, -256, 255, 4, 128, 80, 40); 551 bad_cases (mpfr_compound17, inv_compound17, "mpfr_compound17", 552 0, -256, 255, 4, 128, 80, 40); 553 bad_cases (mpfr_compound120, inv_compound120, "mpfr_compound120", 554 0, -256, 255, 4, 128, 80, 40); 555 556 tests_end_mpfr (); 557 return 0; 558 } 559