1 /* Test file for mpfr_{mul,div}_2{ui,si}.
2
3 Copyright 1999, 2001-2004, 2006-2023 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 char * const val[] = {
26 "1.0001@100","4.0004000000000@102", "4.0004000000000@97",
27 "1.ABF012345@-100","6.afc048d140000@-98","6.afc048d140000@-103",
28 "F.FFFFFFFFF@10000","3.fffffffffc000@10003","3.fffffffffc000@9998",
29 "1.23456789ABCDEF@42","4.8d159e26af37c@44","4.8d159e26af37c@39",
30 "17@42","5.c000000000000@45","5.c000000000000@40",
31 "42@-17","1.0800000000000@-13","1.0800000000000@-18"
32 };
33
34 static int
35 test_mul (int i, int div, mpfr_ptr y, mpfr_srcptr x,
36 unsigned long int n, mpfr_rnd_t r)
37 {
38 return
39 i == 0 ? (div ? mpfr_div_2ui : mpfr_mul_2ui) (y, x, n, r) :
40 i == 1 ? (div ? mpfr_div_2si : mpfr_mul_2si) (y, x, n, r) :
41 i == 2 ? (div ? mpfr_mul_2si : mpfr_div_2si) (y, x, -n, r) :
42 (exit (1), 0);
43 }
44
45 static void
46 underflow (mpfr_exp_t e)
47 {
48 mpfr_t x, y, z1, z2;
49 mpfr_exp_t emin;
50 int i, k, s;
51 int prec;
52 int rnd;
53 int div;
54 int inex1, inex2;
55 unsigned int flags1, flags2;
56
57 /* Test mul_2si(x, e - k), div_2si(x, k - e) and div_2ui(x, k - e) with
58 * emin = e, x = s * (1 + i/16), i in { -1, 0, 1 }, s in { -1, 1 }, and
59 * k = 1 to 4, by comparing the result with the one of a simple division.
60 */
61 emin = mpfr_get_emin ();
62 set_emin (e);
63 mpfr_inits2 (8, x, y, (mpfr_ptr) 0);
64 for (i = 15; i <= 17; i++)
65 for (s = 1; s >= -1; s -= 2)
66 {
67 inex1 = mpfr_set_si_2exp (x, s * i, -4, MPFR_RNDN);
68 MPFR_ASSERTN (inex1 == 0);
69 for (prec = 6; prec >= 3; prec -= 3)
70 {
71 mpfr_inits2 (prec, z1, z2, (mpfr_ptr) 0);
72 RND_LOOP_NO_RNDF (rnd)
73 for (k = 1; k <= 4; k++)
74 {
75 /* The following one is assumed to be correct. */
76 inex1 = mpfr_mul_2si (y, x, e, MPFR_RNDN);
77 MPFR_ASSERTN (inex1 == 0);
78 inex1 = mpfr_set_ui (z1, 1 << k, MPFR_RNDN);
79 MPFR_ASSERTN (inex1 == 0);
80 mpfr_clear_flags ();
81 /* Do not use mpfr_div_ui to avoid the optimization
82 by mpfr_div_2si. */
83 inex1 = mpfr_div (z1, y, z1, (mpfr_rnd_t) rnd);
84 flags1 = __gmpfr_flags;
85
86 for (div = 0; div <= 2; div++)
87 {
88 mpfr_clear_flags ();
89 inex2 =
90 div == 0 ?
91 mpfr_mul_2si (z2, x, e - k, (mpfr_rnd_t) rnd) :
92 div == 1 ?
93 mpfr_div_2si (z2, x, k - e, (mpfr_rnd_t) rnd) :
94 mpfr_div_2ui (z2, x, k - e, (mpfr_rnd_t) rnd);
95 flags2 = __gmpfr_flags;
96 if (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
97 mpfr_equal_p (z1, z2))
98 continue;
99 printf ("Error in underflow(");
100 if (e == MPFR_EMIN_MIN)
101 printf ("MPFR_EMIN_MIN");
102 else if (e == emin)
103 printf ("default emin");
104 else
105 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e);
106 printf (") with mpfr_%s,\nx = %d/16, prec = %d, k = %d,"
107 " %s\n", div == 0 ? "mul_2si" : div == 1 ?
108 "div_2si" : "div_2ui", s * i, prec, k,
109 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
110 printf ("Expected ");
111 mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN);
112 printf (", inex = %d, flags = %u\n",
113 VSIGN (inex1), flags1);
114 printf ("Got ");
115 mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN);
116 printf (", inex = %d, flags = %u\n",
117 VSIGN (inex2), flags2);
118 exit (1);
119 } /* div */
120 } /* k */
121 mpfr_clears (z1, z2, (mpfr_ptr) 0);
122 } /* prec */
123 } /* i */
124 mpfr_clears (x, y, (mpfr_ptr) 0);
125 set_emin (emin);
126 }
127
128 static void
129 underflow0 (void)
130 {
131 underflow (-256);
132 if (mpfr_get_emin () != MPFR_EMIN_MIN)
133 underflow (mpfr_get_emin ());
134 underflow (MPFR_EMIN_MIN);
135 }
136
137 static void
138 large (mpfr_exp_t e)
139 {
140 mpfr_t x, y, z;
141 mpfr_exp_t emax;
142 int inex;
143 unsigned int flags;
144
145 emax = mpfr_get_emax ();
146 set_emax (e);
147 mpfr_init2 (x, 8);
148 mpfr_init2 (y, 8);
149 mpfr_init2 (z, 4);
150
151 mpfr_set_inf (x, 1);
152 mpfr_nextbelow (x);
153
154 mpfr_mul_2si (y, x, -1, MPFR_RNDU);
155 mpfr_prec_round (y, 4, MPFR_RNDU);
156
157 mpfr_clear_flags ();
158 inex = mpfr_mul_2si (z, x, -1, MPFR_RNDU);
159 flags = __gmpfr_flags;
160
161 if (inex <= 0 || flags != MPFR_FLAGS_INEXACT || ! mpfr_equal_p (y, z))
162 {
163 printf ("Error in large(");
164 if (e == MPFR_EMAX_MAX)
165 printf ("MPFR_EMAX_MAX");
166 else if (e == emax)
167 printf ("default emax");
168 else
169 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e);
170 printf (") for mpfr_mul_2si\n");
171 printf ("Expected inex > 0, flags = %u,\n y = ",
172 (unsigned int) MPFR_FLAGS_INEXACT);
173 mpfr_dump (y);
174 printf ("Got inex = %d, flags = %u,\n y = ",
175 inex, flags);
176 mpfr_dump (z);
177 exit (1);
178 }
179
180 mpfr_clear_flags ();
181 inex = mpfr_div_2si (z, x, 1, MPFR_RNDU);
182 flags = __gmpfr_flags;
183
184 if (inex <= 0 || flags != MPFR_FLAGS_INEXACT || ! mpfr_equal_p (y, z))
185 {
186 printf ("Error in large(");
187 if (e == MPFR_EMAX_MAX)
188 printf ("MPFR_EMAX_MAX");
189 else if (e == emax)
190 printf ("default emax");
191 else
192 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e);
193 printf (") for mpfr_div_2si\n");
194 printf ("Expected inex > 0, flags = %u,\n y = ",
195 (unsigned int) MPFR_FLAGS_INEXACT);
196 mpfr_dump (y);
197 printf ("Got inex = %d, flags = %u,\n y = ",
198 inex, flags);
199 mpfr_dump (z);
200 exit (1);
201 }
202
203 mpfr_clear_flags ();
204 inex = mpfr_div_2ui (z, x, 1, MPFR_RNDU);
205 flags = __gmpfr_flags;
206
207 if (inex <= 0 || flags != MPFR_FLAGS_INEXACT || ! mpfr_equal_p (y, z))
208 {
209 printf ("Error in large(");
210 if (e == MPFR_EMAX_MAX)
211 printf ("MPFR_EMAX_MAX");
212 else if (e == emax)
213 printf ("default emax");
214 else
215 printf ("%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e);
216 printf (") for mpfr_div_2ui\n");
217 printf ("Expected inex > 0, flags = %u,\n y = ",
218 (unsigned int) MPFR_FLAGS_INEXACT);
219 mpfr_dump (y);
220 printf ("Got inex = %d, flags = %u,\n y = ",
221 inex, flags);
222 mpfr_dump (z);
223 exit (1);
224 }
225
226 mpfr_clears (x, y, z, (mpfr_ptr) 0);
227 set_emax (emax);
228 }
229
230 static void
231 large0 (void)
232 {
233 mpfr_exp_t emin;
234
235 emin = mpfr_get_emin ();
236
237 while (1)
238 {
239 large (256);
240 if (mpfr_get_emax () != MPFR_EMAX_MAX)
241 large (mpfr_get_emax ());
242 large (MPFR_EMAX_MAX);
243 if (mpfr_get_emin () == MPFR_EMIN_MIN)
244 break;
245 /* Redo the test with __gmpfr_emin set to MPFR_EMIN_MIN, which can
246 be useful to trigger integer overflows as in div_2ui.c r12272. */
247 set_emin (MPFR_EMIN_MIN);
248 }
249
250 set_emin (emin);
251 }
252
253 /* Cases where the function overflows on n = 0 when rounding is like
254 away from zero. */
255 static void
256 overflow0 (mpfr_exp_t emax)
257 {
258 mpfr_exp_t old_emax;
259 mpfr_t x, y1, y2;
260 int neg, r, op;
261 static const char *sop[4] = { "mul_2ui", "mul_2si", "div_2ui", "div_2si" };
262
263 old_emax = mpfr_get_emax ();
264 set_emax (emax);
265
266 mpfr_init2 (x, 8);
267 mpfr_inits2 (6, y1, y2, (mpfr_ptr) 0);
268
269 mpfr_set_inf (x, 1);
270 mpfr_nextbelow (x);
271
272 for (neg = 0; neg <= 1; neg++)
273 {
274 RND_LOOP_NO_RNDF (r)
275 {
276 int inex1, inex2;
277 mpfr_flags_t flags1, flags2;
278
279 /* Even if there isn't an overflow (rounding ~ toward zero),
280 the result is the same as the one of an overflow. */
281 inex1 = mpfr_overflow (y1, (mpfr_rnd_t) r, neg ? -1 : 1);
282 flags1 = MPFR_FLAGS_INEXACT;
283 if (mpfr_inf_p (y1))
284 flags1 |= MPFR_FLAGS_OVERFLOW;
285 for (op = 0; op < 4; op++)
286 {
287 mpfr_clear_flags ();
288 inex2 =
289 op == 0 ? mpfr_mul_2ui (y2, x, 0, (mpfr_rnd_t) r) :
290 op == 1 ? mpfr_mul_2si (y2, x, 0, (mpfr_rnd_t) r) :
291 op == 2 ? mpfr_div_2ui (y2, x, 0, (mpfr_rnd_t) r) :
292 op == 3 ? mpfr_div_2si (y2, x, 0, (mpfr_rnd_t) r) :
293 (MPFR_ASSERTN (0), 0);
294 flags2 = __gmpfr_flags;
295 if (!(mpfr_equal_p (y1, y2) &&
296 SAME_SIGN (inex1, inex2) &&
297 flags1 == flags2))
298 {
299 printf ("Error in overflow0 for %s, mpfr_%s, emax = %"
300 MPFR_EXP_FSPEC "d,\nx = ",
301 mpfr_print_rnd_mode ((mpfr_rnd_t) r), sop[op],
302 (mpfr_eexp_t) emax);
303 mpfr_dump (x);
304 printf ("Expected ");
305 mpfr_dump (y1);
306 printf (" with inex = %d, flags =", inex1);
307 flags_out (flags1);
308 printf ("Got ");
309 mpfr_dump (y2);
310 printf (" with inex = %d, flags =", inex2);
311 flags_out (flags2);
312 exit (1);
313 }
314 }
315 }
316 mpfr_neg (x, x, MPFR_RNDN);
317 }
318
319 mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
320 set_emax (old_emax);
321 }
322
323 static void
324 coverage_div_2ui (void)
325 {
326 mpfr_t x, y;
327
328 mpfr_init2 (x, 2);
329 mpfr_init2 (y, 2);
330 mpfr_set_ui_2exp (x, 1, mpfr_get_emax () - 1, MPFR_RNDN);
331 mpfr_div_2ui (y, x, (unsigned long) LONG_MAX + 1, MPFR_RNDN);
332 MPFR_ASSERTN(mpfr_zero_p (y));
333 MPFR_ASSERTN(mpfr_signbit (y) == 0);
334 mpfr_clear (x);
335 mpfr_clear (y);
336 }
337
338 int
339 main (int argc, char *argv[])
340 {
341 mpfr_t w,z;
342 unsigned long k;
343 int i;
344
345 tests_start_mpfr ();
346
347 coverage_div_2ui ();
348 mpfr_inits2 (53, w, z, (mpfr_ptr) 0);
349
350 for (i = 0; i < 3; i++)
351 {
352 mpfr_set_inf (w, 1);
353 test_mul (i, 0, w, w, 10, MPFR_RNDZ);
354 if (!MPFR_IS_INF(w))
355 {
356 printf ("Result is not Inf (i = %d)\n", i);
357 exit (1);
358 }
359
360 mpfr_set_nan (w);
361 test_mul (i, 0, w, w, 10, MPFR_RNDZ);
362 if (!MPFR_IS_NAN(w))
363 {
364 printf ("Result is not NaN (i = %d)\n", i);
365 exit (1);
366 }
367
368 for (k = 0 ; k < numberof(val) ; k+=3)
369 {
370 mpfr_set_str (w, val[k], 16, MPFR_RNDN);
371 test_mul (i, 0, z, w, 10, MPFR_RNDZ);
372 if (mpfr_cmp_str (z, val[k+1], 16, MPFR_RNDN))
373 {
374 printf ("ERROR for x * 2^n (i = %d) for %s\n", i, val[k]);
375 printf ("Expected: %s\n"
376 "Got : ", val[k+1]);
377 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN);
378 putchar ('\n');
379 exit (1);
380 }
381 test_mul (i, 1, z, w, 10, MPFR_RNDZ);
382 if (mpfr_cmp_str (z, val[k+2], 16, MPFR_RNDN))
383 {
384 printf ("ERROR for x / 2^n (i = %d) for %s\n", i, val[k]);
385 printf ("Expected: %s\n"
386 "Got : ", val[k+2]);
387 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN);
388 putchar ('\n');
389 exit (1);
390 }
391 }
392
393 mpfr_set_inf (w, 1);
394 mpfr_nextbelow (w);
395 test_mul (i, 0, w, w, 1, MPFR_RNDN);
396 if (!mpfr_inf_p (w))
397 {
398 printf ("Overflow error (i = %d)!\n", i);
399 exit (1);
400 }
401 mpfr_set_ui (w, 0, MPFR_RNDN);
402 mpfr_nextabove (w);
403 test_mul (i, 1, w, w, 1, MPFR_RNDN);
404 if (mpfr_cmp_ui (w, 0))
405 {
406 printf ("Underflow error (i = %d)!\n", i);
407 exit (1);
408 }
409 }
410
411 if (MPFR_EXP_MAX >= LONG_MAX/2 && MPFR_EXP_MIN <= LONG_MAX/2-LONG_MAX-1)
412 {
413 unsigned long lmp1 = (unsigned long) LONG_MAX + 1;
414
415 mpfr_set_ui (w, 1, MPFR_RNDN);
416 mpfr_mul_2ui (w, w, LONG_MAX/2, MPFR_RNDZ);
417 mpfr_div_2ui (w, w, lmp1, MPFR_RNDZ);
418 mpfr_mul_2ui (w, w, lmp1 - LONG_MAX/2, MPFR_RNDZ);
419 if (!mpfr_cmp_ui (w, 1))
420 {
421 printf ("Underflow LONG_MAX error!\n");
422 exit (1);
423 }
424 }
425
426 mpfr_clears (w, z, (mpfr_ptr) 0);
427
428 underflow0 ();
429 large0 ();
430
431 if (mpfr_get_emax () != MPFR_EMAX_MAX)
432 overflow0 (mpfr_get_emax ());
433 overflow0 (MPFR_EMAX_MAX);
434 overflow0 (-1);
435
436 tests_end_mpfr ();
437 return 0;
438 }