tests/tsubnormal.c

   1 /* Test file for mpfr_subnormalize.
   2
   3 Copyright 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 #include <stdio.h>
  24 #include <stdlib.h>
  25 #include <limits.h>
  26
  27 #include "mpfr-test.h"
  28
  29 static const struct {
  30   const char *in;
  31   int i;
  32   mpfr_rnd_t rnd;
  33   const char *out;
  34   int j;
  35 } tab[] = { /* 4th field: use the mpfr_dump format, in case of error. */
  36   {"1E1",  0, MPFR_RNDN, "0.100000000E2", 0},
  37   {"1E1", -1, MPFR_RNDZ, "0.100000000E2", -1},
  38   {"1E1", -1, MPFR_RNDD, "0.100000000E2", -1},
  39   {"1E1",  1, MPFR_RNDU, "0.100000000E2", 1},
  40   {"0.10000E-10", 0, MPFR_RNDN, "0.100000000E-10", 0},
  41   {"0.10001E-10", 0, MPFR_RNDN, "0.100000000E-10", -1},
  42   {"0.11001E-10", 0, MPFR_RNDN, "0.100000000E-9", 1},
  43   {"0.11001E-10", 0, MPFR_RNDZ, "0.100000000E-10", -1},
  44   {"0.11001E-10", 0, MPFR_RNDU, "0.100000000E-9", 1},
  45   {"0.11000E-10", 0, MPFR_RNDN, "0.100000000E-9", 1},
  46   {"0.11000E-10", -1, MPFR_RNDN, "0.100000000E-9", 1},
  47   {"0.11000E-10", 1, MPFR_RNDN, "0.100000000E-10", -1},
  48   {"0.11111E-8", 0, MPFR_RNDN, "0.100000000E-7", 1},
  49   {"0.10111E-8", 0, MPFR_RNDN, "0.110000000E-8", 1},
  50   {"0.11110E-8", -1, MPFR_RNDN, "0.100000000E-7", 1},
  51   {"0.10110E-8", 1, MPFR_RNDN, "0.101000000E-8", -1}
  52 };
  53
  54 static void
  55 check1 (void)
  56 {
  57   mpfr_t x;
  58   int i, j, k, s, old_inex, tiny, expj;
  59   mpfr_exp_t emin, emax;
  60   unsigned int expflags, flags;
  61
  62   emin = mpfr_get_emin ();
  63   emax = mpfr_get_emax ();
  64
  65   mpfr_set_default_prec (9);
  66   mpfr_set_emin (-10);
  67   mpfr_set_emax (10);
  68
  69   mpfr_init (x);
  70   for (i = 0; i < (sizeof (tab) / sizeof (tab[0])); i++)
  71     for (s = 0; s <= (tab[i].rnd == MPFR_RNDN); s++)
  72       for (k = 0; k <= 1; k++)
  73         {
  74           mpfr_set_str (x, tab[i].in, 2, MPFR_RNDN);
  75           old_inex = tab[i].i;
  76           expj = tab[i].j;
  77           if (s)
  78             {
  79               mpfr_neg (x, x, MPFR_RNDN);
  80               old_inex = - old_inex;
  81               expj = - expj;
  82             }
  83           if (k && old_inex)
  84             old_inex = old_inex < 0 ? INT_MIN : INT_MAX;
  85           tiny = MPFR_GET_EXP (x) <= -3;
  86           mpfr_clear_flags ();
  87           j = mpfr_subnormalize (x, old_inex, tab[i].rnd);
  88           expflags =
  89             (tiny ? MPFR_FLAGS_UNDERFLOW : 0) |
  90             (expj ? MPFR_FLAGS_INEXACT : 0);
  91           flags = __gmpfr_flags;
  92           if (s)
  93             mpfr_neg (x, x, MPFR_RNDN);
  94           if (mpfr_cmp_str (x, tab[i].out, 2, MPFR_RNDN) != 0 ||
  95               flags != expflags || ! SAME_SIGN (j, expj))
  96             {
  97               const char *sgn = s ? "-" : "";
  98               printf ("Error for i = %d (old_inex = %d), k = %d, x = %s%s\n"
  99                       "Expected: %s%s\nGot:      ", i, old_inex, k,
 100                       sgn, tab[i].in, sgn, tab[i].out);
 101               if (s)
 102                 mpfr_neg (x, x, MPFR_RNDN);
 103               mpfr_dump (x);
 104               printf ("Expected flags = %u, got %u\n", expflags, flags);
 105               printf ("Expected ternary value = %d, got %d\n", expj, j);
 106               exit (1);
 107             }
 108         }
 109   mpfr_clear (x);
 110
 111   MPFR_ASSERTN (mpfr_get_emin () == -10);
 112   MPFR_ASSERTN (mpfr_get_emax () == 10);
 113
 114   set_emin (emin);
 115   set_emax (emax);
 116 }
 117
 118 /* bug found by Kevin P. Rauch on 22 Oct 2007 */
 119 static void
 120 check2 (void)
 121 {
 122   mpfr_t x, y, z;
 123   int tern;
 124   mpfr_exp_t emin;
 125
 126   emin = mpfr_get_emin ();
 127
 128   mpfr_init2 (x, 32);
 129   mpfr_init2 (y, 32);
 130   mpfr_init2 (z, 32);
 131
 132   mpfr_set_ui (x, 0xC0000000U, MPFR_RNDN);
 133   mpfr_neg (x, x, MPFR_RNDN);
 134   mpfr_set_ui (y, 0xFFFFFFFEU, MPFR_RNDN);
 135   mpfr_set_exp (x, 0);
 136   mpfr_set_exp (y, 0);
 137   mpfr_set_emin (-29);
 138
 139   tern = mpfr_mul (z, x, y, MPFR_RNDN);
 140   /* z = -0.BFFFFFFE, tern > 0 */
 141
 142   tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
 143   /* z should be -0.75 */
 144   MPFR_ASSERTN (tern < 0 && mpfr_cmp_si_2exp (z, -3, -2) == 0);
 145
 146   mpfr_clear (x);
 147   mpfr_clear (y);
 148   mpfr_clear (z);
 149
 150   MPFR_ASSERTN (mpfr_get_emin () == -29);
 151
 152   set_emin (emin);
 153 }
 154
 155 /* bug found by Kevin P. Rauch on 22 Oct 2007 */
 156 static void
 157 check3 (void)
 158 {
 159   mpfr_t x, y, z;
 160   int tern;
 161   mpfr_exp_t emin;
 162
 163   emin = mpfr_get_emin ();
 164
 165   mpfr_init2 (x, 32);
 166   mpfr_init2 (y, 32);
 167   mpfr_init2 (z, 32);
 168
 169   mpfr_set_ui (x, 0xBFFFFFFFU, MPFR_RNDN); /* 3221225471/2^32 */
 170   mpfr_set_ui (y, 0x80000001U, MPFR_RNDN); /* 2147483649/2^32 */
 171   mpfr_set_exp (x, 0);
 172   mpfr_set_exp (y, 0);
 173   mpfr_set_emin (-1);
 174
 175   /* the exact product is 6917529028714823679/2^64, which is rounded to
 176      3/8 = 0.375, which is smaller, thus tern < 0 */
 177   tern = mpfr_mul (z, x, y, MPFR_RNDN);
 178   MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0);
 179
 180   tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
 181   /* since emin = -1, and EXP(z)=-1, z should be rounded to precision
 182      EXP(z)-emin+1 = 1, i.e., z should be a multiple of the smallest possible
 183      positive representable value with emin=-1, which is 1/4. The two
 184      possible values are 1/4 and 2/4, which are at equal distance of z.
 185      But since tern < 0, we should choose the largest value, i.e., 2/4. */
 186   MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0);
 187
 188   /* here is another test for the alternate case, where z was rounded up
 189      first, thus we have to round down */
 190   mpfr_set_str_binary (x, "0.11111111111010110101011011011011");
 191   mpfr_set_str_binary (y, "0.01100000000001111100000000001110");
 192   tern = mpfr_mul (z, x, y, MPFR_RNDN);
 193   MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0);
 194   tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
 195   MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 1, -2) == 0);
 196
 197   /* finally the case where z was exact, which we simulate here */
 198   mpfr_set_ui_2exp (z, 3, -3, MPFR_RNDN);
 199   tern = mpfr_subnormalize (z, 0, MPFR_RNDN);
 200   MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0);
 201
 202   mpfr_clear (x);
 203   mpfr_clear (y);
 204   mpfr_clear (z);
 205
 206   MPFR_ASSERTN (mpfr_get_emin () == -1);
 207
 208   set_emin (emin);
 209 }
 210
 211 int
 212 main (int argc, char *argv[])
 213 {
 214   tests_start_mpfr ();
 215
 216   check1 ();
 217   check2 ();
 218   check3 ();
 219
 220   tests_end_mpfr ();
 221   return 0;
 222 }