path: root/tests/tests.c

/* Miscellaneous support for test programs.

Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007 Free Software Foundation, Inc.
Contributed by the Arenaire and Cacao projects, INRIA.

This file is part of the MPFR Library.

The MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the MPFR Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
MA 02110-1301, USA. */

#ifdef HAVE_CONFIG_H
# if HAVE_CONFIG_H
#  include "config.h"     /* for a build within gmp */
# endif
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>

#if HAVE_LOCALE_H
#include <locale.h>
#endif

#if TIME_WITH_SYS_TIME
# include <sys/time.h>  /* for struct timeval */
# include <time.h>
#elif HAVE_SYS_TIME_H
#  include <sys/time.h>
#else
#  include <time.h>
#endif

#if HAVE_SYS_FPU_H
# include <sys/fpu.h>
#endif

#ifdef MPFR_TESTS_TIMEOUT
#include <sys/resource.h>
#endif

#include "mpfr-test.h"

#ifdef MPFR_FPU_PREC
/* This option allows to test MPFR on x86 processors when the FPU
 * rounding precision has been changed. As MPFR is a library, this can
 * occur in practice, either by the calling software or by some other
 * library or plug-in used by the calling software. This option is
 * mainly for developers. If it is used, then the <fpu_control.h>
 * header is assumed to exist and work like under Linux/x86. MPFR does
 * not need to be recompiled. So, a possible usage is the following:
 *
 *   cd tests
 *   make clean
 *   make check CFLAGS="-g -O2 -ffloat-store -DMPFR_FPU_PREC=_FPU_SINGLE"
 *
 * i.e. just add -DMPFR_FPU_PREC=... to the CFLAGS found in Makefile.
 */

#include <fpu_control.h>

static void
set_fpu_prec (void)
{
  fpu_control_t cw;

  _FPU_GETCW(cw);
  cw &= ~(_FPU_EXTENDED|_FPU_DOUBLE|_FPU_SINGLE);
  cw |= (MPFR_FPU_PREC);
  _FPU_SETCW(cw);
}

#endif

static mp_exp_t default_emin, default_emax;

static void tests_rand_start (void);
static void tests_rand_end   (void);
static void tests_limit_start (void);

/* We want to always import the function mpfr_dump inside the test
   suite, so that we can use it in GDB. But it doesn't work if we build
   a Windows DLL (initializer element is not a constant) */
#if !__GMP_LIBGMP_DLL
extern void (*dummy_func) (mpfr_srcptr);
void (*dummy_func)(mpfr_srcptr) = mpfr_dump;
#endif

void
tests_start_mpfr (void)
{
  /* don't buffer, so output is not lost if a test causes a segv etc */
  setbuf (stdout, NULL);

#if HAVE_LOCALE_H && HAVE_SETLOCALE
  /* Added on 2005-07-09. This allows to test MPFR under various
     locales. New bugs will probably be found, in particular with
     LC_ALL="tr_TR.ISO8859-9" because of the i/I character... */
  setlocale (LC_ALL, "");
#endif

#ifdef MPFR_FPU_PREC
  set_fpu_prec ();
#endif

  tests_memory_start ();
  tests_rand_start ();
  tests_limit_start ();

  default_emin = mpfr_get_emin ();
  default_emax = mpfr_get_emax ();
}

void
tests_end_mpfr (void)
{
  int err = 0;

  if (mpfr_get_emin () != default_emin)
    {
      printf ("Default emin value has not been restored!\n");
      err = 1;
    }

  if (mpfr_get_emax () != default_emax)
    {
      printf ("Default emax value has not been restored!\n");
      err = 1;
    }

  mpfr_free_cache ();
  tests_rand_end ();
  tests_memory_end ();

  if (err)
    exit (err);
}

static void
tests_limit_start (void)
{
#ifdef MPFR_TESTS_TIMEOUT
  struct rlimit rlim[1];
  char *timeoutp;
  int timeout;

  timeoutp = getenv ("MPFR_TESTS_TIMEOUT");
  timeout = timeoutp != NULL ? atoi (timeoutp) : MPFR_TESTS_TIMEOUT;
  if (timeout > 0)
    {
      /* We need to call getrlimit first to initialize rlim_max to
         an acceptable value for setrlimit. When enabled, timeouts
         are regarded as important: we don't want to take too much
         CPU time on machines shared with other users. So, if we
         can't set the timeout, we exit immediately. */
      if (getrlimit (RLIMIT_CPU, rlim))
        {
          printf ("Error: getrlimit failed\n");
          exit (1);
        }
      rlim->rlim_cur = timeout;
      if (setrlimit (RLIMIT_CPU, rlim))
        {
          printf ("Error: setrlimit failed\n");
          exit (1);
        }
    }
#endif
}

static void
tests_rand_start (void)
{
  gmp_randstate_ptr  rands;
  char           *perform_seed;
  unsigned long  seed;

  if (__gmp_rands_initialized)
    {
      printf (
        "Please let tests_start() initialize the global __gmp_rands, i.e.\n"
        "ensure that function is called before the first use of RANDS.\n");
      exit (1);
    }

  gmp_randinit_default (__gmp_rands);
  __gmp_rands_initialized = 1;
  rands = __gmp_rands;

  perform_seed = getenv ("GMP_CHECK_RANDOMIZE");
  if (perform_seed != NULL)
    {
      seed = atoi (perform_seed);
      if (! (seed == 0 || seed == 1))
        {
          printf ("Re-seeding with GMP_CHECK_RANDOMIZE=%lu\n", seed);
          gmp_randseed_ui (rands, seed);
        }
      else
        {
#if HAVE_GETTIMEOFDAY
          struct timeval  tv;
          gettimeofday (&tv, NULL);
          seed = tv.tv_sec + tv.tv_usec;
#else
          time_t  tv;
          time (&tv);
          seed = tv;
#endif
          gmp_randseed_ui (rands, seed);
          printf ("Seed GMP_CHECK_RANDOMIZE=%lu "
                  "(include this in bug reports)\n", seed);
        }
    }
  else
      gmp_randseed_ui (rands, 0x2143FEDC);
}

static void
tests_rand_end (void)
{
  RANDS_CLEAR ();
}

/* initialization function for tests using the hardware floats
   Not very useful now. */
void
mpfr_test_init ()
{
  double d;
#if HAVE_FPC_CSR
  /* to get denormalized numbers on IRIX64 */
  union fpc_csr exp;

  exp.fc_word = get_fpc_csr();
  exp.fc_struct.flush = 0;
  set_fpc_csr(exp.fc_word);
#endif
#ifdef HAVE_DENORMS
  d = DBL_MIN;
  if (2.0 * (d / 2.0) != d)
    {
      printf ("Error: HAVE_DENORMS defined, but no subnormals.\n");
      exit (1);
    }
#endif

  /* generate DBL_EPSILON with a loop to avoid that the compiler
     optimizes the code below in non-IEEE 754 mode, deciding that
     c = d is always false. */
#if 0
  for (eps = 1.0; eps != DBL_EPSILON; eps /= 2.0);
  c = 1.0 + eps;
  d = eps * (1.0 - eps) / 2.0;
  d += c;
  if (c != d)
    {
      printf ("Warning: IEEE 754 standard not fully supported\n"
              "         (maybe extended precision not disabled)\n"
              "         Some tests may fail\n");
    }
#endif
}


/* generate a random limb */
mp_limb_t
randlimb (void)
{
  mp_limb_t limb;

  _gmp_rand (&limb, RANDS, BITS_PER_MP_LIMB);
  return limb;
}

/* returns ulp(x) for x a 'normal' double-precision number */
double
Ulp (double x)
{
   double y, eps;

   if (x < 0) x = -x;

   y = x * 2.220446049250313080847263336181640625e-16 ; /* x / 2^52 */

   /* as ulp(x) <= y = x/2^52 < 2*ulp(x),
   we have x + ulp(x) <= x + y <= x + 2*ulp(x),
   therefore o(x + y) = x + ulp(x) or x + 2*ulp(x) */

   eps =  x + y;
   eps = eps - x; /* ulp(x) or 2*ulp(x) */

   return (eps > y) ? 0.5 * eps : eps;
}

/* returns the number of ulp's between a and b,
   where a and b can be any floating-point number, except NaN
 */
int
ulp (double a, double b)
{
  double twoa;

  if (a == b) return 0; /* also deals with a=b=inf or -inf */

  twoa = a + a;
  if (twoa == a) /* a is +/-0.0 or +/-Inf */
    return ((b < a) ? INT_MAX : -INT_MAX);

  return (int) ((a - b) / Ulp (a));
}

/* return double m*2^e */
double
dbl (double m, int e)
{
  if (e >=0 )
    while (e-- > 0)
      m *= 2.0;
  else
    while (e++ < 0)
      m /= 2.0;
  return m;
}

int
Isnan (double d)
{
  return (d) != (d);
}

void
d_trace (const char *name, double d)
{
  union {
    double         d;
    unsigned char  b[sizeof(double)];
  } u;
  int  i;

  if (name != NULL && name[0] != '\0')
    printf ("%s=", name);

  u.d = d;
  printf ("[");
  for (i = 0; i < (int) sizeof (u.b); i++)
    {
      if (i != 0)
        printf (" ");
      printf ("%02X", (int) u.b[i]);
    }
  printf ("] %.20g\n", d);
}

void
ld_trace (const char *name, long double ld)
{
  union {
    long double    ld;
    unsigned char  b[sizeof(long double)];
  } u;
  int  i;

  if (name != NULL && name[0] != '\0')
    printf ("%s=", name);

  u.ld = ld;
  printf ("[");
  for (i = 0; i < (int) sizeof (u.b); i++)
    {
      if (i != 0)
        printf (" ");
      printf ("%02X", (int) u.b[i]);
    }
  printf ("] %.20Lg\n", ld);
}

/* Open a file in the src directory - can't use fopen directly */
FILE *
src_fopen (const char *filename, const char *mode)
{
  const char *srcdir = getenv ("srcdir");
  char *buffer;
  FILE *f;

  if (srcdir == NULL)
    return fopen (filename, mode);
  buffer = (char*) malloc (strlen (filename) + strlen (srcdir) + 2);
  if (buffer == NULL)
    {
      printf ("src_fopen: failed to alloc memory)\n");
      exit (1);
    }
  sprintf (buffer, "%s/%s", srcdir, filename);
  f = fopen (buffer, mode);
  free (buffer);
  return f;
}

void
set_emin (mp_exp_t exponent)
{
  if (mpfr_set_emin (exponent))
    {
      printf ("set_emin: setting emin to %ld failed\n", (long int) exponent);
      exit (1);
    }
}

void
set_emax (mp_exp_t exponent)
{
  if (mpfr_set_emax (exponent))
    {
      printf ("set_emax: setting emax to %ld failed\n", (long int) exponent);
      exit (1);
    }
}

/* pos is 512 times the proportion of negative numbers */
void
tests_default_random (mpfr_ptr x, int pos, mp_exp_t emin, mp_exp_t emax)
{
  mpfr_random (x);
  if (emin >= 1 || (randlimb () & 1))
    mpfr_mul_2si (x, x, emin + (long) (randlimb () % (emax - emin + 1)),
                  GMP_RNDN);
  if (randlimb () % 512 < pos)
    mpfr_neg (x, x, GMP_RNDN);
}

/* The test_one argument is a boolean. If it is true, then the function
   is tested in only one rounding mode (the one provided in rnd) and the
   variable rndnext is not used (due to the break). If it is false, then
   the function is tested in the 4 rounding modes, and rnd must initially
   be GMP_RNDZ; thus rndnext will be initialized in the first iteration.
   gcc may give a warning about rndnext, but this is an easy and correct
   way to implement a simple queue for the rounding modes.
   As examples of use, see the calls to test4rm from the data_check and
   bad_cases functions. */
static void
test4rm (int (*fct) (), mpfr_srcptr x, mpfr_ptr y, mpfr_ptr z,
         mp_rnd_t rnd, int test_one, char *name)
{
  mp_prec_t yprec = MPFR_PREC (y);
  mp_rnd_t rndnext = GMP_RND_MAX;  /* means uninitialized */

  MPFR_ASSERTN (test_one || rnd == GMP_RNDZ);
  mpfr_set_prec (z, yprec);
  while (1)
    {
      MPFR_ASSERTN (rnd != GMP_RND_MAX);
      fct (z, x, rnd);
      if (! mpfr_equal_p (y, z))
        {
          printf ("Error for %s with xprec=%ld, yprec=%ld, rnd=%s\nx = ",
                  name, MPFR_PREC (x), yprec, mpfr_print_rnd_mode (rnd));
          mpfr_out_str (stdout, 16, 0, x, GMP_RNDN);
          printf ("\nexpected ");
          mpfr_out_str (stdout, 16, 0, y, GMP_RNDN);
          printf ("\ngot      ");
          mpfr_out_str (stdout, 16, 0, z, GMP_RNDN);
          printf ("\n");
          exit (1);
        }
      if (test_one || rnd == GMP_RNDN)
        break;
      if (rnd == GMP_RNDZ)
        {
          if (MPFR_IS_NEG (y))
            {
              rnd = GMP_RNDU;
              rndnext = GMP_RNDD;
            }
          else
            {
              rnd = GMP_RNDD;
              rndnext = GMP_RNDU;
            }
        }
      else
        {
          rnd = rndnext;
          if (rndnext != GMP_RNDN)
            {
              rndnext = GMP_RNDN;
              mpfr_nexttoinf (y);
            }
          else
            {
              if (yprec == MPFR_PREC_MIN)
                break;
              mpfr_prec_round (y, --yprec, GMP_RNDZ);
              mpfr_set_prec (z, yprec);
            }
        }
    }
}

/* Check data in file f for function foo, with name 'name'.
   Each line consists of the file f one:

   xprec yprec rnd x y

   where:

   xprec is the input precision
   yprec is the output precision
   rnd is the rounding mode (n, z, u, d, Z)
   x is the input (hexadecimal format)
   y is the expected output (hexadecimal format) for foo(x) with rounding rnd

   If rnd is Z, y is the expected output in round-towards-zero, and the
   three directed rounding modes are tested, then the round-to-nearest
   mode is tested in precision yprec-1. This is useful for worst cases,
   where yprec is the minimum value such that one has a worst case in a
   directed rounding mode.
 */
void
data_check (char *f, int (*foo) (), char *name)
{
  FILE *fp;
  mp_prec_t xprec, yprec;
  mpfr_t x, y, z;
  mp_rnd_t rnd;
  char r;
  int c;

  fp = fopen (f, "r");
  if (fp == NULL)
    {
      char *v = MPFR_VERSION_STRING;

      /* In the '-dev' versions, assume that the data file exists and
         return an error if the file cannot be opened to make sure
         that such failures are detected. */
      while (*v != '\0')
        v++;
      if (v[-4] == '-' && v[-3] == 'd' && v[-2] == 'e' && v[-1] == 'v')
        {
          printf ("Error: unable to open file '%s'\n", f);
          exit (1);
        }
      else
        return;
    }

  mpfr_init (x);
  mpfr_init (y);
  mpfr_init (z);

  while (!feof (fp))
    {
      fscanf (fp, " ");  /* skip whitespace, for consistency */
      c = getc (fp);
      if (c == EOF)
        break;
      if (c == '#') /* comment: read entire line */
        {
          do
            {
              c = getc (fp);
            }
          while (!feof (fp) && c != '\n');
        }
      else
        {
          ungetc (c, fp);
          if (fscanf (fp, "%lu %lu %c", &xprec, &yprec, &r) != 3)
            {
              printf ("Error: corrupted line in file '%s'\n", f);
              exit (1);
            }
          switch (r)
            {
            case 'n':
              rnd = GMP_RNDN;
              break;
            case 'z': case 'Z':
              rnd = GMP_RNDZ;
              break;
            case 'u':
              rnd = GMP_RNDU;
              break;
            case 'd':
              rnd = GMP_RNDD;
              break;
            default:
              printf ("Error: unexpected rounding mode"
                      " in file '%s': %c\n", f, (int) r);
              exit (1);
            }
          mpfr_set_prec (x, xprec);
          mpfr_set_prec (y, yprec);
          if (mpfr_inp_str (x, fp, 0, GMP_RNDN) == 0)
            {
              printf ("Error: corrupted argument in file '%s'\n", f);
              exit (1);
            }
          if (mpfr_inp_str (y, fp, 0, GMP_RNDN) == 0)
            {
              printf ("Error: corrupted result in file '%s'\n", f);
              exit (1);
            }
          if (getc (fp) != '\n')
            {
              printf ("Error: result not followed by \\n in file '%s'\n", f);
              exit (1);
            }
          /* Skip whitespace, in particular at the end of the file. */
          fscanf (fp, " ");
          test4rm (foo, x, y, z, rnd, r != 'Z', name);
        }
    }

  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (z);

  fclose (fp);
}

/* Test n random bad cases. A precision py in [pymin,pymax] and
 * a number y of precision py are chosen randomly. One computes
 * x = inv(y) in precision px = py + psup (rounded to nearest).
 * Then (in general), y is a bad case for fct in precision py (in
 * the directed rounding modes, but also in the rounding-to-nearest
 * mode for some lower precision: see data_check).
 * fct, inv, name: data related to the function.
 * pos, emin, emax: arguments for tests_default_random.
 */
void
bad_cases (int (*fct)(), int (*inv)(), char *name,
           int pos, mp_exp_t emin, mp_exp_t emax,
           mp_prec_t pymin, mp_prec_t pymax, mp_prec_t psup,
           int n)
{
  mpfr_t x, y, z;
  char *dbgenv;
  int i, dbg;

  dbgenv = getenv ("MPFR_DEBUG_BADCASES");
  dbg = dbgenv != 0 ? atoi (dbgenv) : 0;  /* debug level */
  mpfr_inits (x, y, z, (void *) 0);
  for (i = 0; i < n; i++)
    {
      mp_prec_t px, py, pz;
      int inex;

      if (dbg)
        printf ("bad_cases: i = %d\n", i);
      py = pymin + (randlimb () % (pymax - pymin + 1));
      mpfr_set_prec (y, py);
      tests_default_random (y, pos, emin, emax);
      if (dbg)
        {
          printf ("bad_cases: yprec =%4ld, y = ", (long) py);
          mpfr_out_str (stdout, 16, 0, y, GMP_RNDN);
          printf ("\n");
        }
      px = py + psup;
      mpfr_set_prec (x, px);
      mpfr_clear_flags ();
      inv (x, y, GMP_RNDN);
      if (mpfr_nanflag_p () || mpfr_overflow_p () || mpfr_underflow_p ())
        {
          if (dbg)
            printf ("bad_cases: no normal inverse\n");
          goto next_i;
        }
      if (dbg > 1)
        {
          printf ("bad_cases: x = ");
          mpfr_out_str (stdout, 16, 0, x, GMP_RNDN);
          printf ("\n");
        }
      pz = px;
      do
        {
          pz += 32;
          mpfr_set_prec (z, pz);
          if (fct (z, x, GMP_RNDN) == 0)
            {
              if (dbg)
                printf ("bad_cases: exact case\n");
              goto next_i;
            }
          if (dbg)
            {
              if (dbg > 1)
                {
                  printf ("bad_cases: %s(x) ~= ", name);
                  mpfr_out_str (stdout, 16, 0, z, GMP_RNDN);
                }
              else
                {
                  printf ("bad_cases:   [GMP_RNDZ]  ~= ");
                  mpfr_out_str (stdout, 16, 40, z, GMP_RNDZ);
                }
              printf ("\n");
            }
          inex = mpfr_prec_round (z, py, GMP_RNDN);
          if (mpfr_nanflag_p () || mpfr_overflow_p () || mpfr_underflow_p ()
              || ! mpfr_equal_p (z, y))
            {
              if (dbg)
                printf ("bad_cases: inverse doesn't match\n");
              goto next_i;
            }
        }
      while (inex == 0);
      /* We really have a bad case. */
      do
        py--;
      while (py >= MPFR_PREC_MIN && mpfr_prec_round (z, py, GMP_RNDZ) == 0);
      py++;
      /* py is now the smallest output precision such that we have
         a bad case in the directed rounding modes. */
      if (mpfr_prec_round (y, py, GMP_RNDZ) != 0)
        {
          printf ("Internal error for i = %d\n", i);
          exit (1);
        }
      if ((inex > 0 && MPFR_IS_POS (z)) ||
          (inex < 0 && MPFR_IS_NEG (z)))
        {
          mpfr_nexttozero (y);
          if (mpfr_zero_p (y))
            goto next_i;
        }
      if (dbg)
        {
          printf ("bad_cases: yprec =%4ld, y = ", (long) py);
          mpfr_out_str (stdout, 16, 0, y, GMP_RNDN);
          printf ("\n");
        }
      /* Note: y is now the expected result rounded towards zero. */
      test4rm (fct, x, y, z, GMP_RNDZ, 0, name);
    next_i: ;
    }
  mpfr_clears (x, y, z, (void *) 0);
}