tests/tests.c

   1 /* Miscellaneous support for test programs.
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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 http://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 #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
  88 static void
  89 set_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
 101 static mpfr_exp_t default_emin, default_emax;
 102
 103 static void tests_rand_start (void);
 104 static void tests_rand_end   (void);
 105 static 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
 111 extern void (*dummy_func) (mpfr_srcptr);
 112 void (*dummy_func)(mpfr_srcptr) = mpfr_dump;
 113 #endif
 114
 115 void
 116 test_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
 156 void
 157 tests_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
 188 void
 189 tests_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
 227 static void
 228 tests_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
 259 static void
 260 tests_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
 307 static void
 308 tests_rand_end (void)
 309 {
 310   RANDS_CLEAR ();
 311 }
 312
 313 /* initialization function for tests using the hardware floats
 314    Not very useful now. */
 315 void
 316 mpfr_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 */
 355 mp_limb_t
 356 randlimb (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 */
 365 double
 366 Ulp (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  */
 387 int
 388 ulp (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 */
 402 double
 403 dbl (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. */
 415 int
 416 Isnan (double d)
 417 {
 418   return (d) != (d);
 419 }
 420
 421 void
 422 d_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
 444 void
 445 ld_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 */
 468 FILE *
 469 src_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
 493 void
 494 set_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
 503 void
 504 set_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 */
 517 void
 518 tests_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. */
 564 static void
 565 test5rm (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  */
 662 void
 663 data_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  */
 820 void
 821 bad_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
 940 void
 941 flags_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 }