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