1/* Test file for mpfr_subnormalize.
2
3Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by the Arenaire and Cacao 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#include <stdio.h>
24#include <stdlib.h>
25#include <limits.h>
26
27#include "mpfr-test.h"
28
29static 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
54static void
55check1 (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 */
119static void
120check2 (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 */
156static void
157check3 (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
211int
212main (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}
223