tset_ld.c revision 1.1.1.4
1/* Test file for mpfr_set_ld and mpfr_get_ld. 2 3Copyright 2002-2018 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 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include <float.h> 24#ifdef WITH_FPU_CONTROL 25#include <fpu_control.h> 26#endif 27 28#include "mpfr-test.h" 29 30static void 31check_gcc33_bug (void) 32{ 33 volatile long double x; 34 35 x = (long double) 9007199254740992.0 + 1.0; 36 if (x != 0.0) 37 return; /* OK */ 38 printf 39 ("Detected optimization bug of gcc 3.3 on Alpha concerning long double\n" 40 "comparisons; set_ld tests might fail (set_ld won't work correctly).\n" 41 "See https://gcc.gnu.org/ml/gcc-bugs/2003-10/msg00853.html for more\n" 42 "information.\n"); 43} 44 45static int 46Isnan_ld (long double d) 47{ 48 /* Do not convert d to double as this can give an overflow, which 49 may confuse compilers without IEEE 754 support (such as clang 50 -fsanitize=undefined), or trigger a trap if enabled. 51 The DOUBLE_ISNAN macro should work fine on long double. */ 52 if (DOUBLE_ISNAN (d)) 53 return 1; 54 LONGDOUBLE_NAN_ACTION (d, goto yes); 55 return 0; 56 yes: 57 return 1; 58} 59 60/* Return the minimal number of bits to represent d exactly (0 for zero). 61 If flag is non-zero, also print d. */ 62/* FIXME: This function doesn't work if the rounding precision is reduced. */ 63static mpfr_prec_t 64print_binary (long double d, int flag) 65{ 66 long double e, f, r; 67 long exp = 1; 68 mpfr_prec_t prec = 0; 69 70 if (Isnan_ld (d)) 71 { 72 if (flag) 73 printf ("NaN\n"); 74 return 0; 75 } 76 if (d < (long double) 0.0 77#if !defined(MPFR_ERRDIVZERO) 78 || (d == (long double) 0.0 && (1.0 / (double) d < 0.0)) 79#endif 80 ) 81 { 82 if (flag) 83 printf ("-"); 84 d = -d; 85 } 86 /* now d >= 0 */ 87 /* Use 2 differents tests for Inf, to avoid potential bugs 88 in implementations. */ 89 if (Isnan_ld (d - d) || (d > 1 && d * 0.5 == d)) 90 { 91 if (flag) 92 printf ("Inf\n"); 93 return 0; 94 } 95 if (d == (long double) 0.0) 96 { 97 if (flag) 98 printf ("0.0\n"); 99 return prec; 100 } 101 MPFR_ASSERTN (d > 0); 102 e = (long double) 1.0; 103 while (e > d) 104 { 105 e = e * (long double) 0.5; 106 exp --; 107 } 108 if (flag == 2) printf ("1: e=%.36Le\n", e); 109 MPFR_ASSERTN (d >= e); 110 /* FIXME: There can be an overflow here, which may not be supported 111 on all platforms. */ 112 while (f = e + e, d >= f) 113 { 114 e = f; 115 exp ++; 116 } 117 if (flag == 2) printf ("2: e=%.36Le\n", e); 118 MPFR_ASSERTN (e <= d && d < f); 119 if (flag == 1) 120 printf ("0."); 121 if (flag == 2) printf ("3: d=%.36Le e=%.36Le prec=%ld\n", d, e, 122 (long) prec); 123 /* Note: the method we use here to extract the bits of d is the following, 124 to deal with the case where the rounding precision is less than the 125 precision of d: 126 (1) we accumulate the upper bits of d into f 127 (2) when accumulating a new bit into f is not exact, we subtract 128 f from d and reset f to 0 129 This is guaranteed to work only when the rounding precision is at least 130 half the precision of d, since otherwise d-f might not be exact. 131 This method does not work with flush-to-zero on underflow. */ 132 f = 0.0; /* will hold accumulated powers of 2 */ 133 while (1) 134 { 135 prec++; 136 r = f + e; 137 /* r is close to f (in particular in the cases where f+e may 138 not be exact), so that r - f should be exact. */ 139 if (r - f != e) /* f+e is not exact */ 140 { 141 d -= f; /* should be exact */ 142 f = 0.0; 143 r = e; 144 } 145 if (d >= r) 146 { 147 if (flag == 1) 148 printf ("1"); 149 if (d == r) 150 break; 151 f = r; 152 } 153 else 154 { 155 if (flag == 1) 156 printf ("0"); 157 } 158 e *= (long double) 0.5; 159 MPFR_ASSERTN (e != 0); /* may fail with flush-to-zero on underflow */ 160 if (flag == 2) printf ("4: d=%.36Le e=%.36Le prec=%ld\n", d, e, 161 (long) prec); 162 } 163 if (flag == 1) 164 printf ("e%ld\n", exp); 165 return prec; 166} 167 168/* Checks that a long double converted exactly to a MPFR number, then 169 converted back to a long double gives the initial value, or in other 170 words, mpfr_get_ld(mpfr_set_ld(d)) = d. 171*/ 172static void 173check_set_get (long double d) 174{ 175 mpfr_exp_t emin, emax; 176 mpfr_t x; 177 mpfr_prec_t prec; 178 int r; 179 long double e; 180 int inex; 181 int red; 182 183 emin = mpfr_get_emin (); 184 emax = mpfr_get_emax (); 185 186 /* Select a precision to ensure that the conversion of d to x be exact. */ 187 prec = print_binary (d, 0); 188 if (prec < MPFR_PREC_MIN) 189 prec = MPFR_PREC_MIN; 190 mpfr_init2 (x, prec); 191 192 RND_LOOP(r) 193 { 194 inex = mpfr_set_ld (x, d, (mpfr_rnd_t) r); 195 if (inex != 0) 196 { 197 printf ("Error: mpfr_set_ld should be exact (rnd = %s)\n", 198 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 199 /* We use 36 digits here, as the maximum LDBL_MANT_DIG value 200 seen in the current implementations is 113 (binary128), 201 and ceil(1+113*log(2)/log(10)) = 36. But the current glibc 202 implementation of printf with double-double arithmetic 203 (e.g. on PowerPC) is not accurate. */ 204 printf (" d ~= %.36Le (output may be wrong!)\n", d); 205 printf (" inex = %d\n", inex); 206 if (emin >= LONG_MIN) 207 printf (" emin = %ld\n", (long) emin); 208 if (emax <= LONG_MAX) 209 printf (" emax = %ld\n", (long) emax); 210 ld_trace (" d", d); 211 printf (" d = "); 212 print_binary (d, 1); 213 printf (" x = "); 214 mpfr_dump (x); 215 printf (" MPFR_LDBL_MANT_DIG=%u\n", MPFR_LDBL_MANT_DIG); 216 printf (" prec=%lu\n", prec); 217 print_binary (d, 2); 218 exit (1); 219 } 220 for (red = 0; red < 2; red++) 221 { 222 if (red) 223 { 224 mpfr_exp_t ex; 225 226 if (MPFR_IS_SINGULAR (x)) 227 break; 228 ex = MPFR_GET_EXP (x); 229 set_emin (ex); 230 set_emax (ex); 231 } 232 e = mpfr_get_ld (x, (mpfr_rnd_t) r); 233 set_emin (emin); 234 set_emax (emax); 235 if (inex == 0 && ((Isnan_ld(d) && ! Isnan_ld(e)) || 236 (Isnan_ld(e) && ! Isnan_ld(d)) || 237 (e != d && !(Isnan_ld(e) && Isnan_ld(d))))) 238 { 239 printf ("Error: mpfr_get_ld o mpfr_set_ld <> Id%s\n", 240 red ? ", reduced exponent range" : ""); 241 printf (" rnd = %s\n", mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 242 printf (" d ~= %.36Le (output may be wrong!)\n", d); 243 printf (" e ~= %.36Le (output may be wrong!)\n", e); 244 ld_trace (" d", d); 245 printf (" x = "); mpfr_out_str (NULL, 16, 0, x, MPFR_RNDN); 246 printf ("\n"); 247 ld_trace (" e", e); 248 printf (" d = "); 249 print_binary (d, 1); 250 printf (" x = "); 251 mpfr_dump (x); 252 printf (" e = "); 253 print_binary (e, 1); 254 printf (" MPFR_LDBL_MANT_DIG=%u\n", MPFR_LDBL_MANT_DIG); 255#ifdef MPFR_NANISNAN 256 if (Isnan_ld(d) || Isnan_ld(e)) 257 printf ("The reason is that NAN == NAN. Please look at the " 258 "configure output\nand Section \"In case of problem\"" 259 " of the INSTALL file.\n"); 260#endif 261 exit (1); 262 } 263 } 264 } 265 266 mpfr_clear (x); 267} 268 269static void 270test_small (void) 271{ 272 mpfr_t x, y, z; 273 long double d; 274 275 mpfr_init2 (x, MPFR_LDBL_MANT_DIG); 276 mpfr_init2 (y, MPFR_LDBL_MANT_DIG); 277 mpfr_init2 (z, MPFR_LDBL_MANT_DIG); 278 279 /* x = 11906603631607553907/2^(16381+64) */ 280 mpfr_set_str (x, "0.1010010100111100110000001110101101000111010110000001111101110011E-16381", 2, MPFR_RNDN); 281 d = mpfr_get_ld (x, MPFR_RNDN); /* infinite loop? */ 282 mpfr_set_ld (y, d, MPFR_RNDN); 283 mpfr_sub (z, x, y, MPFR_RNDN); 284 mpfr_abs (z, z, MPFR_RNDN); 285 mpfr_clear_erangeflag (); 286 /* If long double = double, d should be equal to 0; 287 in this case, everything is OK. */ 288 if (d != 0 && (mpfr_cmp_str (z, "1E-16434", 2, MPFR_RNDN) > 0 || 289 mpfr_erangeflag_p ())) 290 { 291 printf ("Error with x = "); 292 mpfr_out_str (NULL, 10, 21, x, MPFR_RNDN); 293 printf (" = "); 294 mpfr_out_str (NULL, 16, 0, x, MPFR_RNDN); 295 printf ("\n -> d = %.33Le", d); 296 printf ("\n -> y = "); 297 mpfr_out_str (NULL, 10, 21, y, MPFR_RNDN); 298 printf (" = "); 299 mpfr_out_str (NULL, 16, 0, y, MPFR_RNDN); 300 printf ("\n -> |x-y| = "); 301 mpfr_out_str (NULL, 16, 0, z, MPFR_RNDN); 302 printf ("\n"); 303 exit (1); 304 } 305 306 mpfr_clear (x); 307 mpfr_clear (y); 308 mpfr_clear (z); 309} 310 311static void 312test_fixed_bugs (void) 313{ 314 mpfr_t x; 315 long double l, m; 316 317 /* bug found by Steve Kargl (2009-03-14) */ 318 mpfr_init2 (x, MPFR_LDBL_MANT_DIG); 319 mpfr_set_ui_2exp (x, 1, -16447, MPFR_RNDN); 320 mpfr_get_ld (x, MPFR_RNDN); /* an assertion failed in init2.c:50 */ 321 322 /* bug reported by Jakub Jelinek (2010-10-17) 323 https://gforge.inria.fr/tracker/?func=detail&aid=11300 */ 324 mpfr_set_prec (x, MPFR_LDBL_MANT_DIG); 325 /* l = 0x1.23456789abcdef0123456789abcdp-914L; */ 326 l = 8.215640181713713164092636634579e-276; 327 mpfr_set_ld (x, l, MPFR_RNDN); 328 m = mpfr_get_ld (x, MPFR_RNDN); 329 if (m != l) 330 { 331 printf ("Error in get_ld o set_ld for l=%Le\n", l); 332 printf ("Got m=%Le instead of l\n", m); 333 exit (1); 334 } 335 336 /* another similar test which failed with extended double precision and the 337 generic code for mpfr_set_ld */ 338 /* l = 0x1.23456789abcdef0123456789abcdp-968L; */ 339 l = 4.560596445887084662336528403703e-292; 340 mpfr_set_ld (x, l, MPFR_RNDN); 341 m = mpfr_get_ld (x, MPFR_RNDN); 342 if (m != l) 343 { 344 printf ("Error in get_ld o set_ld for l=%Le\n", l); 345 printf ("Got m=%Le instead of l\n", m); 346 exit (1); 347 } 348 349 mpfr_clear (x); 350} 351 352static void 353check_subnormal (void) 354{ 355 long double d, e; 356 mpfr_t x; 357 358 d = 17.0; 359 mpfr_init2 (x, MPFR_LDBL_MANT_DIG); 360 while (d != 0.0) 361 { 362 mpfr_set_ld (x, d, MPFR_RNDN); 363 e = mpfr_get_ld (x, MPFR_RNDN); 364 if (e != d) 365 { 366 printf ("Error for mpfr_get_ld o mpfr_set_ld\n"); 367 printf ("d=%Le\n", d); 368 printf ("x="); mpfr_dump (x); 369 printf ("e=%Le\n", e); 370 exit (1); 371 } 372 d *= 0.5; 373 } 374 mpfr_clear (x); 375} 376 377static void 378check_overflow (void) 379{ 380 long double d, e; 381 mpfr_t x; 382 int i; 383 384 mpfr_init2 (x, MPFR_LDBL_MANT_DIG); 385 for (i = 0; i < 2; i++) 386 { 387 d = i == 0 ? LDBL_MAX : -LDBL_MAX; 388 mpfr_set_ld (x, d, MPFR_RNDN); 389 mpfr_mul_2ui (x, x, 1, MPFR_RNDN); 390 e = mpfr_get_ld (x, MPFR_RNDN); 391 if (! DOUBLE_ISINF (e) || (i == 0 ? (e < 0) : (e > 0))) 392 { 393 printf ("Error in check_overflow.\n"); 394 printf ("d=%Le\n", d); 395 printf ("x="); mpfr_dump (x); 396 printf ("e=%Le\n", e); 397 exit (1); 398 } 399 } 400 mpfr_clear (x); 401} 402 403/* issue reported by Sisyphus on powerpc */ 404static void 405test_20140212 (void) 406{ 407 mpfr_t fr1, fr2; 408 long double ld, h, l, ld2; 409 int i, c1, c2; 410 411 mpfr_init2 (fr1, 106); 412 mpfr_init2 (fr2, 2098); 413 414 for (h = 1.0L, i = 0; i < 1023; i++) 415 h *= 2.0L; 416 for (l = 1.0L, i = 0; i < 1074; i++) 417 l *= 0.5L; 418 ld = h + l; /* rounding of 2^1023 + 2^(-1074) */ 419 420 mpfr_set_ld (fr1, ld, MPFR_RNDN); 421 mpfr_set_ld (fr2, ld, MPFR_RNDN); 422 423 c1 = mpfr_cmp_ld (fr1, ld); 424 c2 = mpfr_cmp_ld (fr2, ld); 425 426 /* If long double is binary64, then ld = fr1 = fr2 = 2^1023. 427 If long double is double-double, then ld = 2^1023 + 2^(-1074), 428 fr1 = 2^1023 and fr2 = 2^1023 + 2^(-1074) */ 429 MPFR_ASSERTN(ld == h ? (c1 == 0) : (c1 < 0)); 430 431 MPFR_ASSERTN(c2 == 0); 432 433 ld2 = mpfr_get_ld (fr2, MPFR_RNDN); 434 MPFR_ASSERTN(ld2 == ld); 435 436 mpfr_clear (fr1); 437 mpfr_clear (fr2); 438} 439 440/* bug reported by Walter Mascarenhas 441 https://sympa.inria.fr/sympa/arc/mpfr/2016-09/msg00005.html */ 442static void 443bug_20160907 (void) 444{ 445#if HAVE_LDOUBLE_IEEE_EXT_LITTLE 446 long double dn, ld; 447 mpfr_t mp; 448 long e; 449 mpfr_long_double_t x; 450 451 /* the following is the encoding of the smallest subnormal number 452 for HAVE_LDOUBLE_IEEE_EXT_LITTLE */ 453 x.s.manl = 1; 454 x.s.manh = 0; 455 x.s.expl = 0; 456 x.s.exph = 0; 457 x.s.sign= 0; 458 dn = x.ld; 459 e = -16445; 460 /* dn=2^e is now the smallest subnormal. */ 461 462 mpfr_init2 (mp, 64); 463 mpfr_set_ui_2exp (mp, 1, e - 1, MPFR_RNDN); 464 ld = mpfr_get_ld (mp, MPFR_RNDU); 465 /* since mp = 2^(e-1) and ld is rounded upwards, we should have 466 ld = 2^e */ 467 if (ld != dn) 468 { 469 printf ("Error, ld = %Le <> dn = %Le\n", ld, dn); 470 printf ("mp="); 471 mpfr_out_str (stdout, 10, 0, mp, MPFR_RNDN); 472 printf ("\n"); 473 printf ("mp="); mpfr_dump (mp); 474 exit (1); 475 } 476 477 /* check a few more numbers */ 478 for (e = -16446; e <= -16381; e++) 479 { 480 mpfr_set_ui_2exp (mp, 1, e, MPFR_RNDN); 481 ld = mpfr_get_ld (mp, MPFR_RNDU); 482 mpfr_set_ld (mp, ld, MPFR_RNDU); 483 /* mp is 2^e rounded up, thus should be >= 2^e */ 484 MPFR_ASSERTN(mpfr_cmp_ui_2exp (mp, 1, e) >= 0); 485 486 mpfr_set_ui_2exp (mp, 1, e, MPFR_RNDN); 487 ld = mpfr_get_ld (mp, MPFR_RNDD); 488 mpfr_set_ld (mp, ld, MPFR_RNDD); 489 /* mp is 2^e rounded down, thus should be <= 2^e */ 490 if (mpfr_cmp_ui_2exp (mp, 3, e) > 0) 491 { 492 printf ("Error, expected value <= 2^%ld\n", e); 493 printf ("got "); mpfr_dump (mp); 494 exit (1); 495 } 496 } 497 498 mpfr_clear (mp); 499#endif 500} 501 502int 503main (int argc, char *argv[]) 504{ 505 volatile long double d, e, maxp2; 506 mpfr_t x; 507 int i; 508 mpfr_exp_t emax; 509 510 tests_start_mpfr (); 511 mpfr_test_init (); 512 513 check_gcc33_bug (); 514 test_fixed_bugs (); 515 516 mpfr_init2 (x, MPFR_LDBL_MANT_DIG + 64); 517 518#if !defined(MPFR_ERRDIVZERO) 519 /* check NaN */ 520 mpfr_set_nan (x); 521 d = mpfr_get_ld (x, MPFR_RNDN); 522 check_set_get (d); 523#endif 524 525 /* check +0.0 and -0.0 */ 526 d = 0.0; 527 check_set_get (d); 528 d = DBL_NEG_ZERO; 529 check_set_get (d); 530 531 /* check that the sign of -0.0 is set */ 532 mpfr_set_ld (x, DBL_NEG_ZERO, MPFR_RNDN); 533 if (MPFR_IS_POS (x)) 534 { 535#if defined(HAVE_SIGNEDZ) 536 printf ("Error: sign of -0.0 is not set correctly\n"); 537 exit (1); 538#else 539 /* Non IEEE doesn't support negative zero yet */ 540 printf ("Warning: sign of -0.0 is not set correctly\n"); 541#endif 542 } 543 544#if !defined(MPFR_ERRDIVZERO) 545 /* check +Inf */ 546 mpfr_set_inf (x, 1); 547 d = mpfr_get_ld (x, MPFR_RNDN); 548 check_set_get (d); 549 550 /* check -Inf */ 551 mpfr_set_inf (x, -1); 552 d = mpfr_get_ld (x, MPFR_RNDN); 553 check_set_get (d); 554#endif 555 556 /* check the largest power of two */ 557 maxp2 = 1.0; 558 while (maxp2 < LDBL_MAX / 2.0) 559 maxp2 *= 2.0; 560 check_set_get (maxp2); 561 check_set_get (-maxp2); 562 563 d = LDBL_MAX; 564 e = d / 2.0; 565 if (e != maxp2) /* false under NetBSD/x86 */ 566 { 567 /* d = LDBL_MAX does not have excess precision. */ 568 check_set_get (d); 569 check_set_get (-d); 570 } 571 572 /* check the smallest power of two */ 573 d = 1.0; 574 while ((e = d / 2.0) != (long double) 0.0 && e != d) 575 d = e; 576 check_set_get (d); 577 check_set_get (-d); 578 579 /* check that 2^i, 2^i+1, 2^i-1 and 2^i-2^(i-2)-1 are correctly converted */ 580 d = 1.0; 581 for (i = 1; i < MPFR_LDBL_MANT_DIG + 8; i++) 582 { 583 d = 2.0 * d; /* d = 2^i */ 584 check_set_get (d); 585 if (d + 1.0 != d) 586 check_set_get (d + 1.0); 587 else 588 { 589 mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN); 590 mpfr_add_ui (x, x, 1, MPFR_RNDN); 591 e = mpfr_get_ld (x, MPFR_RNDN); 592 check_set_get (e); 593 } 594 if (d - 1.0 != d) 595 check_set_get (d - 1.0); 596 else 597 { 598 mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN); 599 mpfr_sub_ui (x, x, 1, MPFR_RNDN); 600 e = mpfr_get_ld (x, MPFR_RNDN); 601 check_set_get (e); 602 } 603 if (i < 3) 604 continue; 605 /* The following test triggers a failure in r10844 for i = 56, 606 with gcc -mpc64 on x86 (64-bit ABI). */ 607 mpfr_set_ui_2exp (x, 3, i-2, MPFR_RNDN); 608 mpfr_sub_ui (x, x, 1, MPFR_RNDN); 609 e = mpfr_get_ld (x, MPFR_RNDN); 610 check_set_get (e); 611 } 612 613 for (i = 0; i < 10000; i++) 614 { 615 mpfr_urandomb (x, RANDS); 616 d = mpfr_get_ld (x, MPFR_RNDN); 617 check_set_get (d); 618 } 619 620 /* check with reduced emax to exercise overflow */ 621 emax = mpfr_get_emax (); 622 mpfr_set_prec (x, 2); 623 set_emax (1); 624 mpfr_set_ld (x, (long double) 2.0, MPFR_RNDN); 625 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0); 626 for (d = (long double) 2.0, i = 0; i < 13; i++, d *= d); 627 /* now d = 2^8192, or an infinity (e.g. with double or double-double) */ 628 mpfr_set_ld (x, d, MPFR_RNDN); 629 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0); 630 set_emax (emax); 631 632 mpfr_clear (x); 633 634 test_small (); 635 636 check_subnormal (); 637#if !defined(MPFR_ERRDIVZERO) 638 check_overflow (); 639#endif 640 641 test_20140212 (); 642 bug_20160907 (); 643 644 tests_end_mpfr (); 645 646 return 0; 647} 648