1/* Test file for mpfr_subnormalize.
2
3Copyright 2005-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba 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
20https://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 "mpfr-test.h"
24
25static const struct {
26  const char *in;
27  int i;
28  mpfr_rnd_t rnd;
29  const char *out;
30  int j;
31} tab[] = { /* 4th field: use the mpfr_dump format, in case of error. */
32  {"1E1",  0, MPFR_RNDN, "0.100000000E2", 0},
33  {"1E1", -1, MPFR_RNDZ, "0.100000000E2", -1},
34  {"1E1", -1, MPFR_RNDD, "0.100000000E2", -1},
35  {"1E1",  1, MPFR_RNDU, "0.100000000E2", 1},
36  {"0.10000E-10", 0, MPFR_RNDN, "0.100000000E-10", 0},
37  {"0.10001E-10", 0, MPFR_RNDN, "0.100000000E-10", -1},
38  {"0.11001E-10", 0, MPFR_RNDN, "0.100000000E-9", 1},
39  {"0.11001E-10", 0, MPFR_RNDZ, "0.100000000E-10", -1},
40  {"0.11001E-10", 0, MPFR_RNDU, "0.100000000E-9", 1},
41  {"0.11000E-10", 0, MPFR_RNDN, "0.100000000E-9", 1},
42  {"0.11000E-10", -1, MPFR_RNDN, "0.100000000E-9", 1},
43  {"0.11000E-10", 1, MPFR_RNDN, "0.100000000E-10", -1},
44  {"0.11111E-8", 0, MPFR_RNDN, "0.100000000E-7", 1},
45  {"0.10111E-8", 0, MPFR_RNDN, "0.110000000E-8", 1},
46  {"0.11110E-8", -1, MPFR_RNDN, "0.100000000E-7", 1},
47  {"0.10110E-8", 1, MPFR_RNDN, "0.101000000E-8", -1}
48};
49
50static void
51check1 (void)
52{
53  mpfr_t x;
54  int i, j, k, s, old_inex, tiny, expj;
55  mpfr_exp_t emin, emax;
56  unsigned int expflags, flags;
57
58  emin = mpfr_get_emin ();
59  emax = mpfr_get_emax ();
60
61  mpfr_set_default_prec (9);
62  set_emin (-10);
63  set_emax (10);
64
65  mpfr_init (x);
66  for (i = 0; i < numberof (tab); i++)
67    for (s = 0; s <= (tab[i].rnd == MPFR_RNDN); s++)
68      for (k = 0; k <= 1; k++)
69        {
70          mpfr_set_str (x, tab[i].in, 2, MPFR_RNDN);
71          old_inex = tab[i].i;
72          expj = tab[i].j;
73          if (s)
74            {
75              mpfr_neg (x, x, MPFR_RNDN);
76              old_inex = - old_inex;
77              expj = - expj;
78            }
79          if (k && old_inex)
80            old_inex = old_inex < 0 ? INT_MIN : INT_MAX;
81          tiny = MPFR_GET_EXP (x) <= -3;
82          mpfr_clear_flags ();
83          j = mpfr_subnormalize (x, old_inex, tab[i].rnd);
84          expflags =
85            (tiny ? MPFR_FLAGS_UNDERFLOW : 0) |
86            (expj ? MPFR_FLAGS_INEXACT : 0);
87          flags = __gmpfr_flags;
88          if (s)
89            mpfr_neg (x, x, MPFR_RNDN);
90          if (mpfr_cmp_str (x, tab[i].out, 2, MPFR_RNDN) != 0 ||
91              flags != expflags || ! SAME_SIGN (j, expj))
92            {
93              const char *sgn = s ? "-" : "";
94              printf ("Error for i = %d (old_inex = %d), k = %d, x = %s%s\n"
95                      "Expected: %s%s\nGot:      ", i, old_inex, k,
96                      sgn, tab[i].in, sgn, tab[i].out);
97              if (s)
98                mpfr_neg (x, x, MPFR_RNDN);
99              mpfr_dump (x);
100              printf ("Expected flags = %u, got %u\n", expflags, flags);
101              printf ("Expected ternary value = %d, got %d\n", expj, j);
102              exit (1);
103            }
104        }
105  mpfr_clear (x);
106
107  MPFR_ASSERTN (mpfr_get_emin () == -10);
108  MPFR_ASSERTN (mpfr_get_emax () == 10);
109
110  set_emin (emin);
111  set_emax (emax);
112}
113
114/* bug found by Kevin P. Rauch on 22 Oct 2007 */
115static void
116check2 (void)
117{
118  mpfr_t x, y, z;
119  int tern;
120  mpfr_exp_t emin;
121
122  emin = mpfr_get_emin ();
123
124  mpfr_init2 (x, 32);
125  mpfr_init2 (y, 32);
126  mpfr_init2 (z, 32);
127
128  mpfr_set_ui (x, 0xC0000000U, MPFR_RNDN);
129  mpfr_neg (x, x, MPFR_RNDN);
130  mpfr_set_ui (y, 0xFFFFFFFEU, MPFR_RNDN);
131  mpfr_set_exp (x, 0);
132  mpfr_set_exp (y, 0);
133  set_emin (-29);
134
135  tern = mpfr_mul (z, x, y, MPFR_RNDN);
136  /* z = -0.BFFFFFFE, tern > 0 */
137
138  tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
139  /* z should be -0.75 */
140  MPFR_ASSERTN (tern < 0 && mpfr_cmp_si_2exp (z, -3, -2) == 0);
141
142  mpfr_clear (x);
143  mpfr_clear (y);
144  mpfr_clear (z);
145
146  MPFR_ASSERTN (mpfr_get_emin () == -29);
147
148  set_emin (emin);
149}
150
151/* bug found by Kevin P. Rauch on 22 Oct 2007 */
152static void
153check3 (void)
154{
155  mpfr_t x, y, z;
156  int tern;
157  mpfr_exp_t emin;
158
159  emin = mpfr_get_emin ();
160
161  mpfr_init2 (x, 32);
162  mpfr_init2 (y, 32);
163  mpfr_init2 (z, 32);
164
165  mpfr_set_ui (x, 0xBFFFFFFFU, MPFR_RNDN); /* 3221225471/2^32 */
166  mpfr_set_ui (y, 0x80000001U, MPFR_RNDN); /* 2147483649/2^32 */
167  mpfr_set_exp (x, 0);
168  mpfr_set_exp (y, 0);
169  set_emin (-1);
170
171  /* the exact product is 6917529028714823679/2^64, which is rounded to
172     3/8 = 0.375, which is smaller, thus tern < 0 */
173  tern = mpfr_mul (z, x, y, MPFR_RNDN);
174  MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0);
175
176  tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
177  /* since emin = -1, and EXP(z)=-1, z should be rounded to precision
178     EXP(z)-emin+1 = 1, i.e., z should be a multiple of the smallest possible
179     positive representable value with emin=-1, which is 1/4. The two
180     possible values are 1/4 and 2/4, which are at equal distance of z.
181     But since tern < 0, we should choose the largest value, i.e., 2/4. */
182  MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0);
183
184  /* here is another test for the alternate case, where z was rounded up
185     first, thus we have to round down */
186  mpfr_set_str_binary (x, "0.11111111111010110101011011011011");
187  mpfr_set_str_binary (y, "0.01100000000001111100000000001110");
188  tern = mpfr_mul (z, x, y, MPFR_RNDN);
189  MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0);
190  tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
191  MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 1, -2) == 0);
192
193  /* finally the case where z was exact, which we simulate here */
194  mpfr_set_ui_2exp (z, 3, -3, MPFR_RNDN);
195  tern = mpfr_subnormalize (z, 0, MPFR_RNDN);
196  MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0);
197
198  mpfr_clear (x);
199  mpfr_clear (y);
200  mpfr_clear (z);
201
202  MPFR_ASSERTN (mpfr_get_emin () == -1);
203
204  set_emin (emin);
205}
206
207/* check mpfr_subnormize with RNDNA (experimental) */
208static void
209check_rndna (void)
210{
211  mpfr_t x, y;
212  int inex;
213  mpfr_exp_t emin = mpfr_get_emin ();
214
215  mpfr_init2 (x, 11);
216  mpfr_init2 (y, 9);
217
218  mpfr_set_str_binary (x, "0.1111111010E-14");
219  inex = mpfr_set (y, x, MPFR_RNDN);
220  MPFR_ASSERTN(inex == 0);
221  set_emin (-21);
222  inex = mpfr_subnormalize (y, inex, MPFR_RNDNA);
223  /* mpfr_subnormalize rounds to precision EXP(y)-emin+1 = 8,
224     thus should round to 0.111111110E-14 */
225  mpfr_set_str_binary (x, "0.111111110E-14");
226  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
227  MPFR_ASSERTN(inex > 0);
228
229  /* now consider x slightly larger: we should get the same result */
230  mpfr_set_str_binary (x, "0.1111111011E-14");
231  inex = mpfr_set (y, x, MPFR_RNDN);
232  MPFR_ASSERTN(inex > 0);
233  inex = mpfr_subnormalize (y, inex, MPFR_RNDNA);
234  mpfr_set_str_binary (x, "0.111111110E-14");
235  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
236  MPFR_ASSERTN(inex > 0);
237
238  /* now consider x slightly smaller: we should get a different result */
239  mpfr_set_str_binary (x, "0.11111110001E-14");
240  inex = mpfr_set (y, x, MPFR_RNDN);
241  MPFR_ASSERTN(inex < 0);
242  inex = mpfr_subnormalize (y, inex, MPFR_RNDNA);
243  mpfr_set_str_binary (x, "0.111111100E-14");
244  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
245  MPFR_ASSERTN(inex < 0);
246
247  mpfr_clear (x);
248  mpfr_clear (y);
249  set_emin (emin);
250}
251
252/* exercise a corner case of mpfr_subnormalize:
253   y = 1xxx...xxx0|100000| with old_inexact = -1 */
254static void
255coverage (void)
256{
257  mpfr_t y;
258  int inex;
259
260  mpfr_init2 (y, 42);
261  mpfr_set_ui_2exp (y, 131073, mpfr_get_emin () - 2, MPFR_RNDN);
262  MPFR_ASSERTN(mpfr_get_exp (y) == mpfr_get_emin () + 16);
263  /* mpfr_subnormalize rounds y to precision EXP(y) - emin + 1, thus 17 */
264  inex = mpfr_subnormalize (y, -1, MPFR_RNDN);
265  MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 65537, mpfr_get_emin () - 1) == 0);
266  MPFR_ASSERTN(inex > 0);
267  mpfr_clear (y);
268}
269
270int
271main (int argc, char *argv[])
272{
273  tests_start_mpfr ();
274
275  coverage ();
276  check_rndna ();
277  check1 ();
278  check2 ();
279  check3 ();
280
281  tests_end_mpfr ();
282  return 0;
283}
284