1 /* Test file for mpfr_subnormalize.
2 
3 Copyright 2005-2020 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 "mpfr-test.h"
24 
25 static 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 
50 static void
check1(void)51 check1 (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   mpfr_set_emin (-10);
63   mpfr_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 */
115 static void
check2(void)116 check2 (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   mpfr_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 */
152 static void
check3(void)153 check3 (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   mpfr_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) */
208 static void
check_rndna(void)209 check_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   mpfr_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   mpfr_set_emin (emin);
250 }
251 
252 /* exercise a corner case of mpfr_subnormalize:
253    y = 1xxx...xxx0|100000| with old_inexact = -1 */
254 static void
coverage(void)255 coverage (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 
268   mpfr_clear (y);
269 }
270 
271 int
main(int argc,char * argv[])272 main (int argc, char *argv[])
273 {
274   tests_start_mpfr ();
275 
276   coverage ();
277   check_rndna ();
278   check1 ();
279   check2 ();
280   check3 ();
281 
282   tests_end_mpfr ();
283   return 0;
284 }
285