1 /* Test file for mpfr_add1sp.
2
3 Copyright 2004-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 void check_special (void);
26 static void check_random (mpfr_prec_t p);
27
28 static int
29 mpfr_add_cf (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r)
30 {
31 mpfr_clear_flags (); /* allows better checking */
32 return mpfr_add (a, b, c, r);
33 }
34
35 static void
36 check_overflow (void)
37 {
38 mpfr_t x, y, z1, z2;
39 mpfr_exp_t emin, emax;
40
41 emin = mpfr_get_emin ();
42 emax = mpfr_get_emax ();
43
44 set_emin (-1021);
45 set_emax (1024);
46
47 mpfr_inits (x, y, z1, z2, (mpfr_ptr) 0);
48
49 mpfr_set_str1 (x, "8.00468257869324898448e+307");
50 mpfr_set_str1 (y, "7.44784712422708645156e+307");
51 mpfr_add1sp (z1, x, y, MPFR_RNDN);
52 mpfr_add1 (z2, x, y, MPFR_RNDN);
53 if (mpfr_cmp (z1, z2))
54 {
55 printf ("Overflow bug in add1sp.\n");
56 exit (1);
57 }
58 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
59
60 set_emin (emin);
61 set_emax (emax);
62 }
63
64 static void
65 bug20171217 (void)
66 {
67 mpfr_t a, b, c;
68
69 mpfr_init2 (a, 137);
70 mpfr_init2 (b, 137);
71 mpfr_init2 (c, 137);
72 mpfr_set_str_binary (b, "0.11111111111111111111111111111111111111111111111111111111111111111111000000000000000000000000000000000000000000000000000000000000000000000E-66");
73 mpfr_set_str_binary (c, "0.11111111111111111111111111111111111111111111111111111111111111111000000000000000000000000000000000000000000000000000000000000000000110000E-2");
74 mpfr_add_cf (a, b, c, MPFR_RNDN);
75 mpfr_set_str_binary (b, "0.10000000000000000000000000000000000000000000000000000000000000000100000000000000000000000000000000000000000000000000000000000000000001000E-1");
76 MPFR_ASSERTN(mpfr_equal_p (a, b));
77 mpfr_clear (a);
78 mpfr_clear (b);
79 mpfr_clear (c);
80 }
81
82 static void
83 bug20190903 (void)
84 {
85 mpfr_t a, b, c, d;
86 int inex;
87 mpfr_flags_t flags;
88
89 /* Bug in r13574, fixed in r13578.
90 Note: to reproduce the failure, GMP_NUMB_BITS == 64 is assumed. */
91 mpfr_inits2 (128, a, b, c, d, (mpfr_ptr) 0);
92 mpfr_set_str_binary (b, "0.11111111111111111100000000000000000001001111111101111111110000101001111100111110110010011001111110000000101001001001110110101110E0");
93 mpfr_set_str_binary (c, "0.10000001011101000010000111100111011110100000001001000010011011001000110100111111101100001101001101011101100100011000000101110111E-126");
94 mpfr_add_cf (a, b, c, MPFR_RNDN);
95 mpfr_set_str_binary (b, "0.11111111111111111100000000000000000001001111111101111111110000101001111100111110110010011001111110000000101001001001110110110000E0");
96 MPFR_ASSERTN (mpfr_equal_p (a, b));
97 mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
98
99 /* Bug in r13574, fixed in r13586. */
100 /* Figure with GMP_NUMB_BITS = 4:
101 b = 1111 1000
102 c = 1000 0001
103 */
104 mpfr_inits2 (2 * GMP_NUMB_BITS, a, b, c, d, (mpfr_ptr) 0);
105 mpfr_set_ui_2exp (d, 1, 3 * GMP_NUMB_BITS, MPFR_RNDN);
106 mpfr_set_ui_2exp (c, 1, 2 * GMP_NUMB_BITS - 1, MPFR_RNDN);
107 mpfr_sub (b, d, c, MPFR_RNDN);
108 mpfr_add_ui (c, c, 1, MPFR_RNDN);
109 inex = mpfr_add_cf (a, b, c, MPFR_RNDN);
110 flags = __gmpfr_flags;
111 MPFR_ASSERTN (mpfr_equal_p (a, d));
112 MPFR_ASSERTN (inex < 0);
113 MPFR_ASSERTN (flags == MPFR_FLAGS_INEXACT);
114 inex = mpfr_add_cf (a, b, c, MPFR_RNDU);
115 flags = __gmpfr_flags;
116 mpfr_add_ui (d, d, 1, MPFR_RNDU);
117 MPFR_ASSERTN (mpfr_equal_p (a, d));
118 MPFR_ASSERTN (inex > 0);
119 MPFR_ASSERTN (flags == MPFR_FLAGS_INEXACT);
120 mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
121 }
122
123 /* Check corner case b = 1, c = 2^(-p) for MPFR_PREC_MIN <= p <= pmax.
124 With RNDN, result is 1, except for p=1, where it is 2. */
125 static void
126 test_corner_1 (mpfr_prec_t pmax)
127 {
128 mpfr_prec_t p;
129
130 for (p = MPFR_PREC_MIN; p <= pmax; p++)
131 {
132 mpfr_t a, b, c;
133 int inex;
134 mpfr_init2 (a, p);
135 mpfr_init2 (b, p);
136 mpfr_init2 (c, p);
137 mpfr_set_ui (b, 1, MPFR_RNDN);
138 mpfr_set_ui_2exp (c, 1, -p, MPFR_RNDN);
139 inex = mpfr_add_cf (a, b, c, MPFR_RNDN);
140 if (p == 1) /* special case, since 2^(p-1) is odd */
141 {
142 MPFR_ASSERTN(inex > 0);
143 MPFR_ASSERTN(mpfr_cmp_ui (a, 2) == 0);
144 }
145 else
146 {
147 MPFR_ASSERTN(inex < 0);
148 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
149 }
150 mpfr_clear (a);
151 mpfr_clear (b);
152 mpfr_clear (c);
153 }
154 }
155
156 static void
157 coverage (void)
158 {
159 mpfr_t a, b, c;
160 int inex;
161 mpfr_exp_t emax;
162 mpfr_prec_t p;
163
164 mpfr_init (a);
165 mpfr_init (b);
166 mpfr_init (c);
167
168 /* coverage test in mpfr_add1sp: case round away, where add_one_ulp
169 gives a carry, and the new exponent is below emax */
170 for (p = MPFR_PREC_MIN; p <= 3 * GMP_NUMB_BITS; p++)
171 {
172 mpfr_set_prec (a, p);
173 mpfr_set_prec (b, p);
174 mpfr_set_prec (c, p);
175 mpfr_set_ui (b, 1, MPFR_RNDN);
176 mpfr_nextbelow (b); /* b = 1 - 2^(-p) (including for p=1) */
177 mpfr_set_ui_2exp (c, 1, -p-1, MPFR_RNDN);
178 /* c = 2^(-p-1) thus b+c = 1 - 2^(-p-1) should be rounded to 1 */
179 inex = mpfr_add_cf (a, b, c, MPFR_RNDU);
180 MPFR_ASSERTN(inex > 0);
181 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
182 }
183
184 /* coverage test in mpfr_add1sp2: case GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS
185 and a1 = 0 */
186 mpfr_set_prec (a, GMP_NUMB_BITS + 2);
187 mpfr_set_prec (b, GMP_NUMB_BITS + 2);
188 mpfr_set_prec (c, GMP_NUMB_BITS + 2);
189 mpfr_set_ui (b, 1, MPFR_RNDN);
190 mpfr_nextbelow (b); /* b = 1 - 2^(-p) with p = GMP_NUMB_BITS+2 */
191 mpfr_set_ui_2exp (c, 1, -GMP_NUMB_BITS-1, MPFR_RNDN);
192 mpfr_nextbelow (c); /* c = 2^(1-p) - 2^(1-2p) */
193 /* a = 1 + 2^(-p) - 2^(1-2p) should be rounded to 1 with RNDN */
194 inex = mpfr_add_cf (a, b, c, MPFR_RNDN);
195 MPFR_ASSERTN(inex < 0);
196 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
197
198 /* coverage test in mpfr_add1sp2: case round away, where add_one_ulp
199 gives a carry, and the new exponent is below emax */
200 mpfr_set_prec (a, GMP_NUMB_BITS + 1);
201 mpfr_set_prec (b, GMP_NUMB_BITS + 1);
202 mpfr_set_prec (c, GMP_NUMB_BITS + 1);
203 mpfr_set_ui (b, 1, MPFR_RNDN);
204 mpfr_nextbelow (b); /* b = 1 - 2^(-p) */
205 mpfr_set_ui_2exp (c, 1, -GMP_NUMB_BITS-2, MPFR_RNDN);
206 /* c = 2^(-p-1) */
207 inex = mpfr_add_cf (a, b, c, MPFR_RNDU);
208 MPFR_ASSERTN(inex > 0);
209 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
210
211 /* coverage test in mpfr_add1sp3: case GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS
212 and a2 == 0 */
213 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 2);
214 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 2);
215 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 2);
216 mpfr_set_ui (b, 1, MPFR_RNDN);
217 mpfr_nextbelow (b); /* b = 1 - 2^(-p) with p = 2*GMP_NUMB_BITS+2 */
218 mpfr_set_ui_2exp (c, 1, -2*GMP_NUMB_BITS-1, MPFR_RNDN);
219 mpfr_nextbelow (c); /* c = 2^(1-p) - 2^(1-2p) */
220 /* a = 1 + 2^(-p) - 2^(1-2p) should be rounded to 1 with RNDN */
221 inex = mpfr_add_cf (a, b, c, MPFR_RNDN);
222 MPFR_ASSERTN(inex < 0);
223 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
224
225 /* coverage test in mpfr_add1sp3: case bx > emax */
226 emax = mpfr_get_emax ();
227 set_emax (1);
228 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 1);
229 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 1);
230 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 1);
231 mpfr_set_ui_2exp (b, 1, mpfr_get_emax () - 1, MPFR_RNDN);
232 mpfr_nextbelow (b);
233 mpfr_mul_2ui (b, b, 1, MPFR_RNDN);
234 /* now b is the largest number < +Inf */
235 mpfr_div_2ui (c, b, GMP_NUMB_BITS - 1, MPFR_RNDN);
236 /* we are in the case d < GMP_NUMB_BITS of mpfr_add1sp3 */
237 inex = mpfr_add_cf (a, b, b, MPFR_RNDU);
238 MPFR_ASSERTN(inex > 0);
239 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
240 set_emax (emax);
241
242 /* coverage test in mpfr_add1sp3: case round away, where add_one_ulp gives
243 a carry, no overflow */
244 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 1);
245 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 1);
246 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 1);
247 mpfr_set_ui (b, 1, MPFR_RNDN);
248 mpfr_nextbelow (b); /* b = 1 - 2^(-p) */
249 mpfr_set_ui_2exp (c, 1, -2 * GMP_NUMB_BITS - 2, MPFR_RNDN);
250 /* c = 2^(-p-1) */
251 inex = mpfr_add_cf (a, b, c, MPFR_RNDU);
252 MPFR_ASSERTN(inex > 0);
253 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
254
255 /* coverage test in mpfr_add1sp3: case round away, where add_one_ulp gives
256 a carry, with overflow */
257 emax = mpfr_get_emax ();
258 set_emax (1);
259 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 1);
260 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 1);
261 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 1);
262 mpfr_set_ui_2exp (b, 1, mpfr_get_emax () - 1, MPFR_RNDN);
263 mpfr_nextbelow (b);
264 mpfr_mul_2ui (b, b, 1, MPFR_RNDN);
265 /* now b is the largest number < +Inf */
266 mpfr_set_ui_2exp (c, 1, mpfr_get_emin () - 1, MPFR_RNDN);
267 inex = mpfr_add_cf (a, b, c, MPFR_RNDU);
268 MPFR_ASSERTN(inex > 0);
269 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
270 set_emax (emax);
271
272 mpfr_clear (a);
273 mpfr_clear (b);
274 mpfr_clear (c);
275 }
276
277 int
278 main (void)
279 {
280 mpfr_prec_t p;
281 int i;
282
283 tests_start_mpfr ();
284
285 coverage ();
286 test_corner_1 (1024);
287 bug20171217 ();
288 bug20190903 ();
289 check_special ();
290 for (p = MPFR_PREC_MIN; p < 200; p++)
291 check_random (p);
292 for (i = 0; i < 200; i++)
293 {
294 /* special precisions */
295 check_random (GMP_NUMB_BITS);
296 check_random (2 * GMP_NUMB_BITS);
297 }
298 check_overflow ();
299
300 tests_end_mpfr ();
301 return 0;
302 }
303
304 #define STD_ERROR \
305 do \
306 { \
307 printf("ERROR: for %s and p=%lu and i=%d:\nB=", \
308 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \
309 mpfr_dump (b); \
310 printf ("C="); mpfr_dump (c); \
311 printf ("add1 : "); mpfr_dump (a1); \
312 printf ("add1sp: "); mpfr_dump (a2); \
313 exit(1); \
314 } \
315 while (0)
316
317 #define STD_ERROR2 \
318 do \
319 { \
320 printf("ERROR: Wrong inexact flag for %s and p=%lu and i=%d:\nB=", \
321 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \
322 mpfr_dump (b); \
323 printf ("C="); mpfr_dump (c); \
324 printf ("A="); mpfr_dump (a1); \
325 printf ("Add1: %d. Add1sp: %d\n", inexact1, inexact2); \
326 exit(1); \
327 } \
328 while (0)
329
330 #define SET_PREC(_p) \
331 do \
332 { \
333 p = _p; \
334 mpfr_set_prec(a1, _p); mpfr_set_prec(a2, _p); \
335 mpfr_set_prec(b, _p); mpfr_set_prec(c, _p); \
336 } \
337 while (0)
338
339 static void
340 check_random (mpfr_prec_t p)
341 {
342 mpfr_t a1, a2, b, bs, c, cs;
343 int r;
344 int i, inexact1, inexact2;
345
346 mpfr_inits2 (p, a1, a2, b, c, (mpfr_ptr) 0);
347
348 for (i = 0 ; i < 500 ; i++)
349 {
350 mpfr_urandom (b, RANDS, MPFR_RNDA);
351 mpfr_urandom (c, RANDS, MPFR_RNDA);
352 if (MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c))
353 {
354 if (RAND_BOOL ())
355 mpfr_neg (b, b, MPFR_RNDN);
356 if (RAND_BOOL ())
357 mpfr_neg (c, c, MPFR_RNDN);
358 if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c))
359 {
360 /* Exchange b and c, except the signs (actually, the sign
361 of cs doesn't matter). */
362 MPFR_ALIAS (bs, c, MPFR_SIGN (b), MPFR_EXP (c));
363 MPFR_ALIAS (cs, b, MPFR_SIGN (c), MPFR_EXP (b));
364 }
365 else
366 {
367 MPFR_ALIAS (bs, b, MPFR_SIGN (b), MPFR_EXP (b));
368 MPFR_ALIAS (cs, c, MPFR_SIGN (c), MPFR_EXP (c));
369 }
370 RND_LOOP (r)
371 {
372 mpfr_flags_t flags1, flags2;
373
374 if (r == MPFR_RNDF) /* inexact makes no sense, moreover
375 mpfr_add1 and mpfr_add1sp could
376 return different values */
377 continue;
378
379 mpfr_clear_flags ();
380 inexact1 = mpfr_add1 (a1, bs, cs, (mpfr_rnd_t) r);
381 flags1 = __gmpfr_flags;
382 mpfr_clear_flags ();
383 inexact2 = mpfr_add1sp (a2, b, c, (mpfr_rnd_t) r);
384 flags2 = __gmpfr_flags;
385 if (! mpfr_equal_p (a1, a2))
386 STD_ERROR;
387 if (inexact1 != inexact2)
388 STD_ERROR2;
389 MPFR_ASSERTN (flags1 == flags2);
390 }
391 }
392 }
393
394 mpfr_clears (a1, a2, b, c, (mpfr_ptr) 0);
395 }
396
397 static void
398 check_special (void)
399 {
400 mpfr_t a1, a2, b, c;
401 int r;
402 mpfr_prec_t p;
403 int i = -1, inexact1, inexact2;
404
405 mpfr_inits (a1, a2, b, c, (mpfr_ptr) 0);
406
407 RND_LOOP (r)
408 {
409 if (r == MPFR_RNDF)
410 continue; /* inexact makes no sense, mpfr_add1 and mpfr_add1sp
411 could differ */
412 SET_PREC(53);
413 mpfr_set_str1 (b, "1@100");
414 mpfr_set_str1 (c, "1@1");
415 inexact1 = mpfr_add1(a1, b, c, (mpfr_rnd_t) r);
416 inexact2 = mpfr_add1sp(a2, b, c, (mpfr_rnd_t) r);
417 if (mpfr_cmp(a1, a2))
418 STD_ERROR;
419 if (inexact1 != inexact2)
420 STD_ERROR2;
421 mpfr_set_str_binary (b, "1E53");
422 mpfr_set_str_binary (c, "1E0");
423 inexact1 = mpfr_add1(a1, b, c, (mpfr_rnd_t) r);
424 inexact2 = mpfr_add1sp(a2, b, c, (mpfr_rnd_t) r);
425 if (mpfr_cmp(a1, a2))
426 STD_ERROR;
427 if (inexact1 != inexact2)
428 STD_ERROR2;
429 }
430
431 mpfr_set_prec (c, 2);
432 mpfr_set_prec (a1, 2);
433 mpfr_set_prec (a2, 2);
434 mpfr_set_str_binary (c, "1.0e1");
435 mpfr_set_str_binary (a2, "1.1e-1");
436 mpfr_set_str_binary (a1, "0.11E2");
437 mpfr_add1sp (a2, c, a2, MPFR_RNDN);
438 if (mpfr_cmp (a1, a2))
439 {
440 printf ("Regression reuse test failed!\n");
441 exit (1);
442 }
443
444 mpfr_set_prec (a1, 63);
445 mpfr_set_prec (b, 63);
446 mpfr_set_prec (c, 63);
447 mpfr_set_str_binary (b, "0.111111101010110111010100110101010110000101111011011011100011001E-3");
448 mpfr_set_str_binary (c, "0.101111111101110000001100001000011000011011010001010011111100111E-4");
449 mpfr_clear_inexflag ();
450 mpfr_add1sp (a1, b, c, MPFR_RNDN);
451 MPFR_ASSERTN (mpfr_inexflag_p ());
452
453 mpfr_clears (a1, a2, b, c, (mpfr_ptr) 0);
454 }