1/* Miscellaneous support for test programs. 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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#ifdef HAVE_CONFIG_H 24# if HAVE_CONFIG_H 25# include "config.h" /* for a build within gmp */ 26# endif 27#endif 28 29#include <stdlib.h> 30#include <float.h> 31#include <errno.h> 32 33#ifdef HAVE_LOCALE_H 34#include <locale.h> 35#endif 36 37#ifdef TIME_WITH_SYS_TIME 38# include <sys/time.h> /* for struct timeval */ 39# include <time.h> 40#elif defined HAVE_SYS_TIME_H 41# include <sys/time.h> 42#else 43# include <time.h> 44#endif 45 46/* <sys/fpu.h> is needed to have union fpc_csr defined under IRIX64 47 (see below). Let's include it only if need be. */ 48#if defined HAVE_SYS_FPU_H && defined HAVE_FPC_CSR 49# include <sys/fpu.h> 50#endif 51 52#ifdef MPFR_TESTS_TIMEOUT 53#include <sys/resource.h> 54#endif 55 56#include "mpfr-test.h" 57 58#ifdef MPFR_FPU_PREC 59/* This option allows to test MPFR on x86 processors when the FPU 60 * rounding precision has been changed. As MPFR is a library, this can 61 * occur in practice, either by the calling software or by some other 62 * library or plug-in used by the calling software. This option is 63 * mainly for developers. If it is used, then the <fpu_control.h> 64 * header is assumed to exist and work like under Linux/x86. MPFR does 65 * not need to be recompiled. So, a possible usage is the following: 66 * 67 * cd tests 68 * make clean 69 * make check CFLAGS="-g -O2 -ffloat-store -DMPFR_FPU_PREC=_FPU_SINGLE" 70 * 71 * i.e. just add -DMPFR_FPU_PREC=... to the CFLAGS found in Makefile. 72 * 73 * Notes: 74 * + SSE2 (used to implement double's on x86_64, and possibly on x86 75 * too, depending on the compiler configuration and flags) is not 76 * affected by the dynamic precision. 77 * + When the FPU is set to single precision, the behavior of MPFR 78 * functions that have a native floating-point type (float, double, 79 * long double) as argument or return value is not guaranteed. 80 */ 81 82#include <fpu_control.h> 83 84static void 85set_fpu_prec (void) 86{ 87 fpu_control_t cw; 88 89 _FPU_GETCW(cw); 90 cw &= ~(_FPU_EXTENDED|_FPU_DOUBLE|_FPU_SINGLE); 91 cw |= (MPFR_FPU_PREC); 92 _FPU_SETCW(cw); 93} 94 95#endif 96 97static mpfr_exp_t default_emin, default_emax; 98 99static void tests_rand_start (void); 100static void tests_rand_end (void); 101static void tests_limit_start (void); 102 103/* We want to always import the function mpfr_dump inside the test 104 suite, so that we can use it in GDB. But it doesn't work if we build 105 a Windows DLL (initializer element is not a constant) */ 106#if !__GMP_LIBGMP_DLL 107extern void (*dummy_func) (mpfr_srcptr); 108void (*dummy_func)(mpfr_srcptr) = mpfr_dump; 109#endif 110 111void 112test_version (void) 113{ 114 const char *version; 115 116 /* VL: I get the following error on an OpenSUSE machine, and changing 117 the value of shlibpath_overrides_runpath in the libtool file from 118 'no' to 'yes' fixes the problem. */ 119 120 version = mpfr_get_version (); 121 if (strcmp (MPFR_VERSION_STRING, version) == 0) 122 { 123 char buffer[16]; 124 int i; 125 126 sprintf (buffer, "%d.%d.%d", MPFR_VERSION_MAJOR, MPFR_VERSION_MINOR, 127 MPFR_VERSION_PATCHLEVEL); 128 for (i = 0; buffer[i] == version[i]; i++) 129 if (buffer[i] == '\0') 130 return; 131 if (buffer[i] == '\0' && version[i] == '-') 132 return; 133 printf ("MPFR_VERSION_MAJOR.MPFR_VERSION_MINOR.MPFR_VERSION_PATCHLEVEL" 134 " (%s)\nand MPFR_VERSION_STRING (%s) do not match!\nIt seems " 135 "that the mpfr.h file has been corrupted.\n", buffer, version); 136 exit (1); 137 } 138 139 printf ("Incorrect MPFR version! (%s header vs %s library)\n" 140 "Nothing else has been tested since for this reason,\n" 141 "any other test may fail. Please fix this one first.\n\n" 142 "You can try to avoid this problem by changing the value of\n" 143 "shlibpath_overrides_runpath in the libtool file and rebuild\n" 144 "MPFR (make clean && make && make check).\n" 145 "Otherwise this error may be due to a corrupted mpfr.h, an\n" 146 "incomplete build (try to rebuild MPFR from scratch and/or\n" 147 "use 'make clean'), or something wrong in the system.\n", 148 MPFR_VERSION_STRING, version); 149 exit (1); 150} 151 152void 153tests_start_mpfr (void) 154{ 155 test_version (); 156 157 /* don't buffer, so output is not lost if a test causes a segv etc */ 158 setbuf (stdout, NULL); 159 160#if defined HAVE_LOCALE_H && defined HAVE_SETLOCALE 161 /* Added on 2005-07-09. This allows to test MPFR under various 162 locales. New bugs will probably be found, in particular with 163 LC_ALL="tr_TR.ISO8859-9" because of the i/I character... */ 164 setlocale (LC_ALL, ""); 165#endif 166 167#ifdef MPFR_FPU_PREC 168 set_fpu_prec (); 169#endif 170 171 tests_memory_start (); 172 tests_rand_start (); 173 tests_limit_start (); 174 175 default_emin = mpfr_get_emin (); 176 default_emax = mpfr_get_emax (); 177} 178 179void 180tests_end_mpfr (void) 181{ 182 int err = 0; 183 184 if (mpfr_get_emin () != default_emin) 185 { 186 printf ("Default emin value has not been restored!\n"); 187 err = 1; 188 } 189 190 if (mpfr_get_emax () != default_emax) 191 { 192 printf ("Default emax value has not been restored!\n"); 193 err = 1; 194 } 195 196 mpfr_free_cache (); 197 tests_rand_end (); 198 tests_memory_end (); 199 200 if (err) 201 exit (err); 202} 203 204static void 205tests_limit_start (void) 206{ 207#ifdef MPFR_TESTS_TIMEOUT 208 struct rlimit rlim[1]; 209 char *timeoutp; 210 int timeout; 211 212 timeoutp = getenv ("MPFR_TESTS_TIMEOUT"); 213 timeout = timeoutp != NULL ? atoi (timeoutp) : MPFR_TESTS_TIMEOUT; 214 if (timeout > 0) 215 { 216 /* We need to call getrlimit first to initialize rlim_max to 217 an acceptable value for setrlimit. When enabled, timeouts 218 are regarded as important: we don't want to take too much 219 CPU time on machines shared with other users. So, if we 220 can't set the timeout, we exit immediately. */ 221 if (getrlimit (RLIMIT_CPU, rlim)) 222 { 223 printf ("Error: getrlimit failed\n"); 224 exit (1); 225 } 226 rlim->rlim_cur = timeout; 227 if (setrlimit (RLIMIT_CPU, rlim)) 228 { 229 printf ("Error: setrlimit failed\n"); 230 exit (1); 231 } 232 } 233#endif 234} 235 236static void 237tests_rand_start (void) 238{ 239 gmp_randstate_ptr rands; 240 char *perform_seed; 241 unsigned long seed; 242 243 if (__gmp_rands_initialized) 244 { 245 printf ( 246 "Please let tests_start() initialize the global __gmp_rands, i.e.\n" 247 "ensure that function is called before the first use of RANDS.\n"); 248 exit (1); 249 } 250 251 gmp_randinit_default (__gmp_rands); 252 __gmp_rands_initialized = 1; 253 rands = __gmp_rands; 254 255 perform_seed = getenv ("GMP_CHECK_RANDOMIZE"); 256 if (perform_seed != NULL) 257 { 258 seed = atoi (perform_seed); 259 if (! (seed == 0 || seed == 1)) 260 { 261 printf ("Re-seeding with GMP_CHECK_RANDOMIZE=%lu\n", seed); 262 gmp_randseed_ui (rands, seed); 263 } 264 else 265 { 266#ifdef HAVE_GETTIMEOFDAY 267 struct timeval tv; 268 gettimeofday (&tv, NULL); 269 seed = tv.tv_sec + tv.tv_usec; 270#else 271 time_t tv; 272 time (&tv); 273 seed = tv; 274#endif 275 gmp_randseed_ui (rands, seed); 276 printf ("Seed GMP_CHECK_RANDOMIZE=%lu " 277 "(include this in bug reports)\n", seed); 278 } 279 } 280 else 281 gmp_randseed_ui (rands, 0x2143FEDC); 282} 283 284static void 285tests_rand_end (void) 286{ 287 RANDS_CLEAR (); 288} 289 290/* initialization function for tests using the hardware floats 291 Not very useful now. */ 292void 293mpfr_test_init (void) 294{ 295 double d; 296#ifdef HAVE_FPC_CSR 297 /* to get denormalized numbers on IRIX64 */ 298 union fpc_csr exp; 299 300 exp.fc_word = get_fpc_csr(); 301 exp.fc_struct.flush = 0; 302 set_fpc_csr(exp.fc_word); 303#endif 304#ifdef HAVE_DENORMS 305 d = DBL_MIN; 306 if (2.0 * (d / 2.0) != d) 307 { 308 printf ("Error: HAVE_DENORMS defined, but no subnormals.\n"); 309 exit (1); 310 } 311#endif 312 313 /* generate DBL_EPSILON with a loop to avoid that the compiler 314 optimizes the code below in non-IEEE 754 mode, deciding that 315 c = d is always false. */ 316#if 0 317 for (eps = 1.0; eps != DBL_EPSILON; eps /= 2.0); 318 c = 1.0 + eps; 319 d = eps * (1.0 - eps) / 2.0; 320 d += c; 321 if (c != d) 322 { 323 printf ("Warning: IEEE 754 standard not fully supported\n" 324 " (maybe extended precision not disabled)\n" 325 " Some tests may fail\n"); 326 } 327#endif 328} 329 330 331/* generate a random limb */ 332mp_limb_t 333randlimb (void) 334{ 335 mp_limb_t limb; 336 337 mpfr_rand_raw (&limb, RANDS, GMP_NUMB_BITS); 338 return limb; 339} 340 341/* returns ulp(x) for x a 'normal' double-precision number */ 342double 343Ulp (double x) 344{ 345 double y, eps; 346 347 if (x < 0) x = -x; 348 349 y = x * 2.220446049250313080847263336181640625e-16 ; /* x / 2^52 */ 350 351 /* as ulp(x) <= y = x/2^52 < 2*ulp(x), 352 we have x + ulp(x) <= x + y <= x + 2*ulp(x), 353 therefore o(x + y) = x + ulp(x) or x + 2*ulp(x) */ 354 355 eps = x + y; 356 eps = eps - x; /* ulp(x) or 2*ulp(x) */ 357 358 return (eps > y) ? 0.5 * eps : eps; 359} 360 361/* returns the number of ulp's between a and b, 362 where a and b can be any floating-point number, except NaN 363 */ 364int 365ulp (double a, double b) 366{ 367 double twoa; 368 369 if (a == b) return 0; /* also deals with a=b=inf or -inf */ 370 371 twoa = a + a; 372 if (twoa == a) /* a is +/-0.0 or +/-Inf */ 373 return ((b < a) ? INT_MAX : -INT_MAX); 374 375 return (int) ((a - b) / Ulp (a)); 376} 377 378/* return double m*2^e */ 379double 380dbl (double m, int e) 381{ 382 if (e >=0 ) 383 while (e-- > 0) 384 m *= 2.0; 385 else 386 while (e++ < 0) 387 m /= 2.0; 388 return m; 389} 390 391/* Warning: NaN values cannot be distinguished if MPFR_NANISNAN is defined. */ 392int 393Isnan (double d) 394{ 395 return (d) != (d); 396} 397 398void 399d_trace (const char *name, double d) 400{ 401 union { 402 double d; 403 unsigned char b[sizeof(double)]; 404 } u; 405 int i; 406 407 if (name != NULL && name[0] != '\0') 408 printf ("%s=", name); 409 410 u.d = d; 411 printf ("["); 412 for (i = 0; i < (int) sizeof (u.b); i++) 413 { 414 if (i != 0) 415 printf (" "); 416 printf ("%02X", (int) u.b[i]); 417 } 418 printf ("] %.20g\n", d); 419} 420 421void 422ld_trace (const char *name, long double ld) 423{ 424 union { 425 long double ld; 426 unsigned char b[sizeof(long double)]; 427 } u; 428 int i; 429 430 if (name != NULL && name[0] != '\0') 431 printf ("%s=", name); 432 433 u.ld = ld; 434 printf ("["); 435 for (i = 0; i < (int) sizeof (u.b); i++) 436 { 437 if (i != 0) 438 printf (" "); 439 printf ("%02X", (int) u.b[i]); 440 } 441 printf ("] %.20Lg\n", ld); 442} 443 444/* Open a file in the src directory - can't use fopen directly */ 445FILE * 446src_fopen (const char *filename, const char *mode) 447{ 448 const char *srcdir = getenv ("srcdir"); 449 char *buffer; 450 FILE *f; 451 452 if (srcdir == NULL) 453 return fopen (filename, mode); 454 buffer = 455 (char*) (*__gmp_allocate_func) (strlen (filename) + strlen (srcdir) + 2); 456 if (buffer == NULL) 457 { 458 printf ("src_fopen: failed to alloc memory)\n"); 459 exit (1); 460 } 461 sprintf (buffer, "%s/%s", srcdir, filename); 462 f = fopen (buffer, mode); 463 (*__gmp_free_func) (buffer, strlen (filename) + strlen (srcdir) + 2); 464 return f; 465} 466 467void 468set_emin (mpfr_exp_t exponent) 469{ 470 if (mpfr_set_emin (exponent)) 471 { 472 printf ("set_emin: setting emin to %ld failed\n", (long int) exponent); 473 exit (1); 474 } 475} 476 477void 478set_emax (mpfr_exp_t exponent) 479{ 480 if (mpfr_set_emax (exponent)) 481 { 482 printf ("set_emax: setting emax to %ld failed\n", (long int) exponent); 483 exit (1); 484 } 485} 486 487/* pos is 512 times the proportion of negative numbers. 488 If pos=256, half of the numbers are negative. 489 If pos=0, all generated numbers are positive. 490*/ 491void 492tests_default_random (mpfr_ptr x, int pos, mpfr_exp_t emin, mpfr_exp_t emax) 493{ 494 MPFR_ASSERTN (emin <= emax); 495 MPFR_ASSERTN (emin >= MPFR_EMIN_MIN); 496 MPFR_ASSERTN (emax <= MPFR_EMAX_MAX); 497 /* but it isn't required that emin and emax are in the current 498 exponent range (see below), so that underflow/overflow checks 499 can be done on 64-bit machines. */ 500 501 mpfr_urandomb (x, RANDS); 502 if (MPFR_IS_PURE_FP (x) && (emin >= 1 || (randlimb () & 1))) 503 { 504 mpfr_exp_t e; 505 e = MPFR_GET_EXP (x) + 506 (emin + (long) (randlimb () % (emax - emin + 1))); 507 /* Note: There should be no overflow here because both terms are 508 between MPFR_EMIN_MIN and MPFR_EMAX_MAX, but the sum e isn't 509 necessarily between MPFR_EMIN_MIN and MPFR_EMAX_MAX. */ 510 if (mpfr_set_exp (x, e)) 511 { 512 /* The random number doesn't fit in the current exponent range. 513 In this case, test the function in the extended exponent range, 514 which should be restored by the caller. */ 515 mpfr_set_emin (MPFR_EMIN_MIN); 516 mpfr_set_emax (MPFR_EMAX_MAX); 517 mpfr_set_exp (x, e); 518 } 519 } 520 if (randlimb () % 512 < pos) 521 mpfr_neg (x, x, MPFR_RNDN); 522} 523 524/* The test_one argument is seen a boolean. If it is true and rnd is 525 a rounding mode toward infinity, then the function is tested in 526 only one rounding mode (the one provided in rnd) and the variable 527 rndnext is not used (due to the break). If it is true and rnd is a 528 rounding mode toward or away from zero, then the function is tested 529 twice, first with the provided rounding mode and second with the 530 rounding mode toward the corresponding infinity (determined by the 531 sign of the result). If it is false, then the function is tested 532 in the 5 rounding modes, and rnd must initially be MPFR_RNDZ; thus 533 rndnext will be initialized in the first iteration. 534 If the test_one argument is 2, then this means that y is exact, and 535 the ternary value is checked. 536 As examples of use, see the calls to test5rm from the data_check and 537 bad_cases functions. */ 538static void 539test5rm (int (*fct) (FLIST), mpfr_srcptr x, mpfr_ptr y, mpfr_ptr z, 540 mpfr_rnd_t rnd, int test_one, char *name) 541{ 542 mpfr_prec_t yprec = MPFR_PREC (y); 543 mpfr_rnd_t rndnext = MPFR_RND_MAX; /* means uninitialized */ 544 545 MPFR_ASSERTN (test_one || rnd == MPFR_RNDZ); 546 mpfr_set_prec (z, yprec); 547 while (1) 548 { 549 int inex; 550 551 MPFR_ASSERTN (rnd != MPFR_RND_MAX); 552 inex = fct (z, x, rnd); 553 if (! (mpfr_equal_p (y, z) || (mpfr_nan_p (y) && mpfr_nan_p (z)))) 554 { 555 printf ("Error for %s with xprec=%lu, yprec=%lu, rnd=%s\nx = ", 556 name, (unsigned long) MPFR_PREC (x), (unsigned long) yprec, 557 mpfr_print_rnd_mode (rnd)); 558 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); 559 printf ("\nexpected "); 560 mpfr_out_str (stdout, 16, 0, y, MPFR_RNDN); 561 printf ("\ngot "); 562 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN); 563 printf ("\n"); 564 exit (1); 565 } 566 if (test_one == 2 && inex != 0) 567 { 568 printf ("Error for %s with xprec=%lu, yprec=%lu, rnd=%s\nx = ", 569 name, (unsigned long) MPFR_PREC (x), (unsigned long) yprec, 570 mpfr_print_rnd_mode (rnd)); 571 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); 572 printf ("\nexact case, but non-zero ternary value (%d)\n", inex); 573 exit (1); 574 } 575 if (rnd == MPFR_RNDN) 576 break; 577 578 if (test_one) 579 { 580 if (rnd == MPFR_RNDU || rnd == MPFR_RNDD) 581 break; 582 583 if (MPFR_IS_NEG (y)) 584 rnd = (rnd == MPFR_RNDA) ? MPFR_RNDD : MPFR_RNDU; 585 else 586 rnd = (rnd == MPFR_RNDA) ? MPFR_RNDU : MPFR_RNDD; 587 } 588 else if (rnd == MPFR_RNDZ) 589 { 590 rnd = MPFR_IS_NEG (y) ? MPFR_RNDU : MPFR_RNDD; 591 rndnext = MPFR_RNDA; 592 } 593 else 594 { 595 rnd = rndnext; 596 if (rnd == MPFR_RNDA) 597 { 598 mpfr_nexttoinf (y); 599 rndnext = (MPFR_IS_NEG (y)) ? MPFR_RNDD : MPFR_RNDU; 600 } 601 else if (rndnext != MPFR_RNDN) 602 rndnext = MPFR_RNDN; 603 else 604 { 605 if (yprec == MPFR_PREC_MIN) 606 break; 607 mpfr_prec_round (y, --yprec, MPFR_RNDZ); 608 mpfr_set_prec (z, yprec); 609 } 610 } 611 } 612} 613 614/* Check data in file f for function foo, with name 'name'. 615 Each line consists of the file f one: 616 617 xprec yprec rnd x y 618 619 where: 620 621 xprec is the input precision 622 yprec is the output precision 623 rnd is the rounding mode (n, z, u, d, a, Z, *) 624 x is the input (hexadecimal format) 625 y is the expected output (hexadecimal format) for foo(x) with rounding rnd 626 627 If rnd is Z, y is the expected output in round-toward-zero, and the 628 four directed rounding modes are tested, then the round-to-nearest 629 mode is tested in precision yprec-1. This is useful for worst cases, 630 where yprec is the minimum value such that one has a worst case in a 631 directed rounding mode. 632 633 If rnd is *, y must be an exact case. All the rounding modes are tested 634 and the ternary value is checked (it must be 0). 635 */ 636void 637data_check (char *f, int (*foo) (FLIST), char *name) 638{ 639 FILE *fp; 640 int xprec, yprec; /* not mpfr_prec_t because of the fscanf */ 641 mpfr_t x, y, z; 642 mpfr_rnd_t rnd; 643 char r; 644 int c; 645 646 fp = fopen (f, "r"); 647 if (fp == NULL) 648 fp = src_fopen (f, "r"); 649 if (fp == NULL) 650 { 651 char *v = (char *) MPFR_VERSION_STRING; 652 653 /* In the '-dev' versions, assume that the data file exists and 654 return an error if the file cannot be opened to make sure 655 that such failures are detected. */ 656 while (*v != '\0') 657 v++; 658 if (v[-4] == '-' && v[-3] == 'd' && v[-2] == 'e' && v[-1] == 'v') 659 { 660 printf ("Error: unable to open file '%s'\n", f); 661 exit (1); 662 } 663 else 664 return; 665 } 666 667 mpfr_init (x); 668 mpfr_init (y); 669 mpfr_init (z); 670 671 while (!feof (fp)) 672 { 673 /* skip whitespace, for consistency */ 674 if (fscanf (fp, " ") == EOF) 675 { 676 if (ferror (fp)) 677 { 678 perror ("data_check"); 679 exit (1); 680 } 681 else 682 break; /* end of file */ 683 } 684 685 if ((c = getc (fp)) == EOF) 686 { 687 if (ferror (fp)) 688 { 689 perror ("data_check"); 690 exit (1); 691 } 692 else 693 break; /* end of file */ 694 } 695 696 if (c == '#') /* comment: read entire line */ 697 { 698 do 699 { 700 c = getc (fp); 701 } 702 while (!feof (fp) && c != '\n'); 703 } 704 else 705 { 706 ungetc (c, fp); 707 708 c = fscanf (fp, "%d %d %c", &xprec, &yprec, &r); 709 MPFR_ASSERTN (xprec >= MPFR_PREC_MIN && xprec <= MPFR_PREC_MAX); 710 MPFR_ASSERTN (yprec >= MPFR_PREC_MIN && yprec <= MPFR_PREC_MAX); 711 if (c == EOF) 712 { 713 perror ("data_check"); 714 exit (1); 715 } 716 else if (c != 3) 717 { 718 printf ("Error: corrupted line in file '%s'\n", f); 719 exit (1); 720 } 721 722 switch (r) 723 { 724 case 'n': 725 rnd = MPFR_RNDN; 726 break; 727 case 'z': case 'Z': 728 rnd = MPFR_RNDZ; 729 break; 730 case 'u': 731 rnd = MPFR_RNDU; 732 break; 733 case 'd': 734 rnd = MPFR_RNDD; 735 break; 736 case '*': 737 rnd = MPFR_RND_MAX; /* non-existing rounding mode */ 738 break; 739 default: 740 printf ("Error: unexpected rounding mode" 741 " in file '%s': %c\n", f, (int) r); 742 exit (1); 743 } 744 mpfr_set_prec (x, xprec); 745 mpfr_set_prec (y, yprec); 746 if (mpfr_inp_str (x, fp, 0, MPFR_RNDN) == 0) 747 { 748 printf ("Error: corrupted argument in file '%s'\n", f); 749 exit (1); 750 } 751 if (mpfr_inp_str (y, fp, 0, MPFR_RNDN) == 0) 752 { 753 printf ("Error: corrupted result in file '%s'\n", f); 754 exit (1); 755 } 756 if (getc (fp) != '\n') 757 { 758 printf ("Error: result not followed by \\n in file '%s'\n", f); 759 exit (1); 760 } 761 /* Skip whitespace, in particular at the end of the file. */ 762 if (fscanf (fp, " ") == EOF && ferror (fp)) 763 { 764 perror ("data_check"); 765 exit (1); 766 } 767 if (r == '*') 768 { 769 int rndint; 770 RND_LOOP (rndint) 771 test5rm (foo, x, y, z, (mpfr_rnd_t) rndint, 2, name); 772 } 773 else 774 test5rm (foo, x, y, z, rnd, r != 'Z', name); 775 } 776 } 777 778 mpfr_clear (x); 779 mpfr_clear (y); 780 mpfr_clear (z); 781 782 fclose (fp); 783} 784 785/* Test n random bad cases. A precision py in [pymin,pymax] and 786 * a number y of precision py are chosen randomly. One computes 787 * x = inv(y) in precision px = py + psup (rounded to nearest). 788 * Then (in general), y is a bad case for fct in precision py (in 789 * the directed rounding modes, but also in the rounding-to-nearest 790 * mode for some lower precision: see data_check). 791 * fct, inv, name: data related to the function. 792 * pos, emin, emax: arguments for tests_default_random. 793 */ 794void 795bad_cases (int (*fct)(FLIST), int (*inv)(FLIST), char *name, 796 int pos, mpfr_exp_t emin, mpfr_exp_t emax, 797 mpfr_prec_t pymin, mpfr_prec_t pymax, mpfr_prec_t psup, 798 int n) 799{ 800 mpfr_t x, y, z; 801 char *dbgenv; 802 int i, dbg; 803 mpfr_exp_t old_emin, old_emax; 804 805 old_emin = mpfr_get_emin (); 806 old_emax = mpfr_get_emax (); 807 808 dbgenv = getenv ("MPFR_DEBUG_BADCASES"); 809 dbg = dbgenv != 0 ? atoi (dbgenv) : 0; /* debug level */ 810 mpfr_inits (x, y, z, (mpfr_ptr) 0); 811 for (i = 0; i < n; i++) 812 { 813 mpfr_prec_t px, py, pz; 814 int inex; 815 816 if (dbg) 817 printf ("bad_cases: i = %d\n", i); 818 py = pymin + (randlimb () % (pymax - pymin + 1)); 819 mpfr_set_prec (y, py); 820 tests_default_random (y, pos, emin, emax); 821 if (dbg) 822 { 823 printf ("bad_cases: yprec =%4ld, y = ", (long) py); 824 mpfr_out_str (stdout, 16, 0, y, MPFR_RNDN); 825 printf ("\n"); 826 } 827 px = py + psup; 828 mpfr_set_prec (x, px); 829 mpfr_clear_flags (); 830 inv (x, y, MPFR_RNDN); 831 if (mpfr_nanflag_p () || mpfr_overflow_p () || mpfr_underflow_p ()) 832 { 833 if (dbg) 834 printf ("bad_cases: no normal inverse\n"); 835 goto next_i; 836 } 837 if (dbg > 1) 838 { 839 printf ("bad_cases: x = "); 840 mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); 841 printf ("\n"); 842 } 843 pz = px; 844 do 845 { 846 pz += 32; 847 mpfr_set_prec (z, pz); 848 if (fct (z, x, MPFR_RNDN) == 0) 849 { 850 if (dbg) 851 printf ("bad_cases: exact case\n"); 852 goto next_i; 853 } 854 if (dbg) 855 { 856 if (dbg > 1) 857 { 858 printf ("bad_cases: %s(x) ~= ", name); 859 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN); 860 } 861 else 862 { 863 printf ("bad_cases: [MPFR_RNDZ] ~= "); 864 mpfr_out_str (stdout, 16, 40, z, MPFR_RNDZ); 865 } 866 printf ("\n"); 867 } 868 inex = mpfr_prec_round (z, py, MPFR_RNDN); 869 if (mpfr_nanflag_p () || mpfr_overflow_p () || mpfr_underflow_p () 870 || ! mpfr_equal_p (z, y)) 871 { 872 if (dbg) 873 printf ("bad_cases: inverse doesn't match\n"); 874 goto next_i; 875 } 876 } 877 while (inex == 0); 878 /* We really have a bad case. */ 879 do 880 py--; 881 while (py >= MPFR_PREC_MIN && mpfr_prec_round (z, py, MPFR_RNDZ) == 0); 882 py++; 883 /* py is now the smallest output precision such that we have 884 a bad case in the directed rounding modes. */ 885 if (mpfr_prec_round (y, py, MPFR_RNDZ) != 0) 886 { 887 printf ("Internal error for i = %d\n", i); 888 exit (1); 889 } 890 if ((inex > 0 && MPFR_IS_POS (z)) || 891 (inex < 0 && MPFR_IS_NEG (z))) 892 { 893 mpfr_nexttozero (y); 894 if (mpfr_zero_p (y)) 895 goto next_i; 896 } 897 if (dbg) 898 { 899 printf ("bad_cases: yprec =%4ld, y = ", (long) py); 900 mpfr_out_str (stdout, 16, 0, y, MPFR_RNDN); 901 printf ("\n"); 902 } 903 /* Note: y is now the expected result rounded toward zero. */ 904 test5rm (fct, x, y, z, MPFR_RNDZ, 0, name); 905 next_i: 906 /* In case the exponent range has been changed by 907 tests_default_random()... */ 908 mpfr_set_emin (old_emin); 909 mpfr_set_emax (old_emax); 910 } 911 mpfr_clears (x, y, z, (mpfr_ptr) 0); 912} 913