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