path: root/tests/trint.c

/* Test file for mpfr_rint, mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round.

Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
Contributed by the Arenaire and Caramel projects, INRIA.

This file is part of the GNU MPFR Library.

The GNU 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 3 of the License, or (at your
option) any later version.

The GNU 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 GNU MPFR Library; see the file COPYING.LESSER.  If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */

#include <stdlib.h>

#include "mpfr-test.h"

#if __MPFR_STDC (199901L)
# include <math.h>
#endif

static void
special (void)
{
  mpfr_t x, y;
  mpfr_exp_t emax;

  mpfr_init (x);
  mpfr_init (y);

  mpfr_set_nan (x);
  mpfr_rint (y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_nan_p (y));

  mpfr_set_inf (x, 1);
  mpfr_rint (y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0);

  mpfr_set_inf (x, -1);
  mpfr_rint (y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) < 0);

  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_rint (y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y));

  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_neg (x, x, MPFR_RNDN);
  mpfr_rint (y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_NEG(y));

  /* coverage test */
  mpfr_set_prec (x, 2);
  mpfr_set_ui (x, 1, MPFR_RNDN);
  mpfr_mul_2exp (x, x, mp_bits_per_limb, MPFR_RNDN);
  mpfr_rint (y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);

  /* another coverage test */
  emax = mpfr_get_emax ();
  set_emax (1);
  mpfr_set_prec (x, 3);
  mpfr_set_str_binary (x, "1.11E0");
  mpfr_set_prec (y, 2);
  mpfr_rint (y, x, MPFR_RNDU); /* x rounds to 1.0E1=0.1E2 which overflows */
  MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0);
  set_emax (emax);

  /* yet another */
  mpfr_set_prec (x, 97);
  mpfr_set_prec (y, 96);
  mpfr_set_str_binary (x, "-0.1011111001101111000111011100011100000110110110110000000111010001000101001111101010101011010111100E97");
  mpfr_rint (y, x, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);

  mpfr_set_prec (x, 53);
  mpfr_set_prec (y, 53);
  mpfr_set_str_binary (x, "0.10101100000000101001010101111111000000011111010000010E-1");
  mpfr_rint (y, x, MPFR_RNDU);
  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);
  mpfr_rint (y, x, MPFR_RNDD);
  MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y));

  mpfr_set_prec (x, 36);
  mpfr_set_prec (y, 2);
  mpfr_set_str_binary (x, "-11000110101010111111110111001.0000100");
  mpfr_rint (y, x, MPFR_RNDN);
  mpfr_set_str_binary (x, "-11E27");
  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);

  mpfr_set_prec (x, 39);
  mpfr_set_prec (y, 29);
  mpfr_set_str_binary (x, "-0.100010110100011010001111001001001100111E39");
  mpfr_rint (y, x, MPFR_RNDN);
  mpfr_set_str_binary (x, "-0.10001011010001101000111100101E39");
  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);

  mpfr_set_prec (x, 46);
  mpfr_set_prec (y, 32);
  mpfr_set_str_binary (x, "-0.1011100110100101000001011111101011001001101001E32");
  mpfr_rint (y, x, MPFR_RNDN);
  mpfr_set_str_binary (x, "-0.10111001101001010000010111111011E32");
  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);

  /* coverage test for mpfr_round */
  mpfr_set_prec (x, 3);
  mpfr_set_str_binary (x, "1.01E1"); /* 2.5 */
  mpfr_set_prec (y, 2);
  mpfr_round (y, x);
  /* since mpfr_round breaks ties away, should give 3 and not 2 as with
     the "round to even" rule */
  MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0);
  /* same test for the function */
  (mpfr_round) (y, x);
  MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0);

  mpfr_set_prec (x, 6);
  mpfr_set_prec (y, 3);
  mpfr_set_str_binary (x, "110.111");
  mpfr_round (y, x);
  if (mpfr_cmp_ui (y, 7))
    {
      printf ("Error in round(110.111)\n");
      exit (1);
    }

  /* Bug found by  Mark J Watkins */
  mpfr_set_prec (x, 84);
  mpfr_set_str_binary (x,
   "0.110011010010001000000111101101001111111100101110010000000000000" \
                       "000000000000000000000E32");
  mpfr_round (x, x);
  if (mpfr_cmp_str (x, "0.1100110100100010000001111011010100000000000000" \
                    "00000000000000000000000000000000000000E32", 2, MPFR_RNDN))
    {
      printf ("Rounding error when dest=src\n");
      exit (1);
    }

  mpfr_clear (x);
  mpfr_clear (y);
}

#if __MPFR_STDC (199901L)

static void
test_fct (double (*f)(double), int (*g)(), char *s, mpfr_rnd_t r)
{
  double d, y;
  mpfr_t dd, yy;

  mpfr_init2 (dd, 53);
  mpfr_init2 (yy, 53);
  for (d = -5.0; d <= 5.0; d += 0.25)
    {
      mpfr_set_d (dd, d, r);
      y = (*f)(d);
      if (g == &mpfr_rint)
        mpfr_rint (yy, dd, r);
      else
        (*g)(yy, dd);
      mpfr_set_d (dd, y, r);
      if (mpfr_cmp (yy, dd))
        {
          printf ("test_against_libc: incorrect result for %s, rnd = %s,"
                  " d = %g\ngot ", s, mpfr_print_rnd_mode (r), d);
          mpfr_out_str (stdout, 10, 0, yy, MPFR_RNDN);
          printf (" instead of %g\n", y);
          exit (1);
        }
    }
  mpfr_clear (dd);
  mpfr_clear (yy);
}

#define TEST_FCT(F) test_fct (&F, &mpfr_##F, #F, r)

static void
test_against_libc (void)
{
  mpfr_rnd_t r = MPFR_RNDN;

#if HAVE_ROUND
  TEST_FCT (round);
#endif
#if HAVE_TRUNC
  TEST_FCT (trunc);
#endif
#if HAVE_FLOOR
  TEST_FCT (floor);
#endif
#if HAVE_CEIL
  TEST_FCT (ceil);
#endif
#if HAVE_NEARBYINT
  for (r = 0; r < MPFR_RND_MAX ; r++)
    if (mpfr_set_machine_rnd_mode (r) == 0)
      test_fct (&nearbyint, &mpfr_rint, "rint", r);
#endif
}

#endif

static void
err (const char *str, mp_size_t s, mpfr_t x, mpfr_t y, mpfr_prec_t p,
     mpfr_rnd_t r, int trint, int inexact)
{
  printf ("Error: %s\ns = %u, p = %u, r = %s, trint = %d, inexact = %d\nx = ",
          str, (unsigned int) s, (unsigned int) p, mpfr_print_rnd_mode (r),
          trint, inexact);
  mpfr_print_binary (x);
  printf ("\ny = ");
  mpfr_print_binary (y);
  printf ("\n");
  exit (1);
}

static void
coverage_03032011 (void)
{
  mpfr_t in, out, cmp;
  int status;
  int precIn;
  char strData[(GMP_NUMB_BITS * 4)+256];

  precIn = GMP_NUMB_BITS * 4;

  mpfr_init2 (in, precIn);
  mpfr_init2 (out, GMP_NUMB_BITS);
  mpfr_init2 (cmp, GMP_NUMB_BITS);

  /* cmp = "0.1EprecIn+2" */
  /* The buffer size is sufficient, as precIn is small in practice. */
  sprintf (strData, "0.1E%d", precIn+2);
  mpfr_set_str_binary (cmp, strData);

  /* in = "0.10...01EprecIn+2" use all (precIn) significand bits */
  memset ((void *)strData, '0', precIn+2);
  strData[1] = '.';
  strData[2] = '1';
  sprintf (&strData[precIn+1], "1E%d", precIn+2);
  mpfr_set_str_binary (in, strData);

  status = mpfr_rint (out, in, MPFR_RNDN);
  if ((mpfr_cmp (out, cmp) != 0) || (status >= 0))
    {
      printf("mpfr_rint error :\n status is %d instead of 0\n", status);
      printf(" out value is ");
      mpfr_dump(out);
      printf(" instead of   ");
      mpfr_dump(cmp);
      exit (1);
    }

  mpfr_clear (cmp);
  mpfr_clear (out);

  mpfr_init2 (out, GMP_NUMB_BITS);
  mpfr_init2 (cmp, GMP_NUMB_BITS);

  /* cmp = "0.10...01EprecIn+2" use all (GMP_NUMB_BITS) significand bits */
  strcpy (&strData[GMP_NUMB_BITS+1], &strData[precIn+1]);
  mpfr_set_str_binary (cmp, strData);

  (MPFR_MANT(in))[2] = MPFR_LIMB_HIGHBIT;
  status = mpfr_rint (out, in, MPFR_RNDN);

  if ((mpfr_cmp (out, cmp) != 0) || (status <= 0))
    {
      printf("mpfr_rint error :\n status is %d instead of 0\n", status);
      printf(" out value is\n");
      mpfr_dump(out);
      printf(" instead of\n");
      mpfr_dump(cmp);
      exit (1);
    }

  mpfr_clear (cmp);
  mpfr_clear (out);
  mpfr_clear (in);
}

int
main (int argc, char *argv[])
{
  mp_size_t s;
  mpz_t z;
  mpfr_prec_t p;
  mpfr_t x, y, t, u, v;
  int r;
  int inexact, sign_t;

  tests_start_mpfr ();

  mpfr_init (x);
  mpfr_init (y);
  mpz_init (z);
  mpfr_init (t);
  mpfr_init (u);
  mpfr_init (v);
  mpz_set_ui (z, 1);
  for (s = 2; s < 100; s++)
    {
      /* z has exactly s bits */

      mpz_mul_2exp (z, z, 1);
      if (randlimb () % 2)
        mpz_add_ui (z, z, 1);
      mpfr_set_prec (x, s);
      mpfr_set_prec (t, s);
      mpfr_set_prec (u, s);
      if (mpfr_set_z (x, z, MPFR_RNDN))
        {
          printf ("Error: mpfr_set_z should be exact (s = %u)\n",
                  (unsigned int) s);
          exit (1);
        }
      if (randlimb () % 2)
        mpfr_neg (x, x, MPFR_RNDN);
      if (randlimb () % 2)
        mpfr_div_2ui (x, x, randlimb () % s, MPFR_RNDN);
      for (p = 2; p < 100; p++)
        {
          int trint;
          mpfr_set_prec (y, p);
          mpfr_set_prec (v, p);
          for (r = 0; r < MPFR_RND_MAX ; r++)
            for (trint = 0; trint < 3; trint++)
              {
                if (trint == 2)
                  inexact = mpfr_rint (y, x, (mpfr_rnd_t) r);
                else if (r == MPFR_RNDN)
                  inexact = mpfr_round (y, x);
                else if (r == MPFR_RNDZ)
                  inexact = (trint ? mpfr_trunc (y, x) :
                             mpfr_rint_trunc (y, x, MPFR_RNDZ));
                else if (r == MPFR_RNDU)
                  inexact = (trint ? mpfr_ceil (y, x) :
                             mpfr_rint_ceil (y, x, MPFR_RNDU));
                else /* r = MPFR_RNDD */
                  inexact = (trint ? mpfr_floor (y, x) :
                             mpfr_rint_floor (y, x, MPFR_RNDD));
                if (mpfr_sub (t, y, x, MPFR_RNDN))
                  err ("subtraction 1 should be exact",
                       s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                sign_t = mpfr_cmp_ui (t, 0);
                if (trint != 0 &&
                    (((inexact == 0) && (sign_t != 0)) ||
                     ((inexact < 0) && (sign_t >= 0)) ||
                     ((inexact > 0) && (sign_t <= 0))))
                  err ("wrong inexact flag", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                if (inexact == 0)
                  continue; /* end of the test for exact results */

                if (((r == MPFR_RNDD || (r == MPFR_RNDZ && MPFR_SIGN (x) > 0))
                     && inexact > 0) ||
                    ((r == MPFR_RNDU || (r == MPFR_RNDZ && MPFR_SIGN (x) < 0))
                     && inexact < 0))
                  err ("wrong rounding direction",
                       s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                if (inexact < 0)
                  {
                    mpfr_add_ui (v, y, 1, MPFR_RNDU);
                    if (mpfr_cmp (v, x) <= 0)
                      err ("representable integer between x and its "
                           "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                  }
                else
                  {
                    mpfr_sub_ui (v, y, 1, MPFR_RNDD);
                    if (mpfr_cmp (v, x) >= 0)
                      err ("representable integer between x and its "
                           "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                  }
                if (r == MPFR_RNDN)
                  {
                    int cmp;
                    if (mpfr_sub (u, v, x, MPFR_RNDN))
                      err ("subtraction 2 should be exact",
                           s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                    cmp = mpfr_cmp_abs (t, u);
                    if (cmp > 0)
                      err ("faithful rounding, but not the nearest integer",
                           s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                    if (cmp < 0)
                      continue;
                    /* |t| = |u|: x is the middle of two consecutive
                       representable integers. */
                    if (trint == 2)
                      {
                        /* halfway case for mpfr_rint in MPFR_RNDN rounding
                           mode: round to an even integer or significand. */
                        mpfr_div_2ui (y, y, 1, MPFR_RNDZ);
                        if (!mpfr_integer_p (y))
                          err ("halfway case for mpfr_rint, result isn't an"
                               " even integer", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                        /* If floor(x) and ceil(x) aren't both representable
                           integers, the significand must be even. */
                        mpfr_sub (v, v, y, MPFR_RNDN);
                        mpfr_abs (v, v, MPFR_RNDN);
                        if (mpfr_cmp_ui (v, 1) != 0)
                          {
                            mpfr_div_2si (y, y, MPFR_EXP (y) - MPFR_PREC (y)
                                          + 1, MPFR_RNDN);
                            if (!mpfr_integer_p (y))
                              err ("halfway case for mpfr_rint, significand isn't"
                                   " even", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                          }
                      }
                    else
                      { /* halfway case for mpfr_round: x must have been
                           rounded away from zero. */
                        if ((MPFR_SIGN (x) > 0 && inexact < 0) ||
                            (MPFR_SIGN (x) < 0 && inexact > 0))
                          err ("halfway case for mpfr_round, bad rounding"
                               " direction", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
                      }
                  }
              }
        }
    }
  mpfr_clear (x);
  mpfr_clear (y);
  mpz_clear (z);
  mpfr_clear (t);
  mpfr_clear (u);
  mpfr_clear (v);

  special ();
  coverage_03032011 ();

#if __MPFR_STDC (199901L)
  if (argc > 1 && strcmp (argv[1], "-s") == 0)
    test_against_libc ();
#endif

  tests_end_mpfr ();
  return 0;
}