1 /* Test file for mpfr_fmma and mpfr_fmms.
2
3 Copyright 2016-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 /* check both mpfr_fmma and mpfr_fmms */
26 static void
27 random_test (mpfr_ptr a, mpfr_ptr b, mpfr_ptr c, mpfr_ptr d, mpfr_rnd_t rnd)
28 {
29 mpfr_t ref, res, ab, cd;
30 int inex_ref, inex_res;
31 mpfr_prec_t p = MPFR_PREC(a);
32
33 mpfr_init2 (res, p);
34 mpfr_init2 (ref, p);
35 mpfr_init2 (ab, mpfr_get_prec (a) + mpfr_get_prec (b));
36 mpfr_init2 (cd, mpfr_get_prec (c) + mpfr_get_prec (d));
37
38 /* first check fmma */
39 inex_res = mpfr_fmma (res, a, b, c, d, rnd);
40 mpfr_mul (ab, a, b, rnd);
41 mpfr_mul (cd, c, d, rnd);
42 inex_ref = mpfr_add (ref, ab, cd, rnd);
43 if (! SAME_SIGN (inex_res, inex_ref) ||
44 mpfr_nan_p (res) || mpfr_nan_p (ref) ||
45 ! mpfr_equal_p (res, ref))
46 {
47 printf ("mpfr_fmma failed for p=%ld rnd=%s\n",
48 (long int) p, mpfr_print_rnd_mode (rnd));
49 printf ("a="); mpfr_dump (a);
50 printf ("b="); mpfr_dump (b);
51 printf ("c="); mpfr_dump (c);
52 printf ("d="); mpfr_dump (d);
53 printf ("Expected\n ");
54 mpfr_dump (ref);
55 printf (" with inex = %d\n", inex_ref);
56 printf ("Got\n ");
57 mpfr_dump (res);
58 printf (" with inex = %d\n", inex_res);
59 exit (1);
60 }
61
62 /* then check fmms */
63 inex_res = mpfr_fmms (res, a, b, c, d, rnd);
64 mpfr_mul (ab, a, b, rnd);
65 mpfr_mul (cd, c, d, rnd);
66 inex_ref = mpfr_sub (ref, ab, cd, rnd);
67 if (! SAME_SIGN (inex_res, inex_ref) ||
68 mpfr_nan_p (res) || mpfr_nan_p (ref) ||
69 ! mpfr_equal_p (res, ref))
70 {
71 printf ("mpfr_fmms failed for p=%ld rnd=%s\n",
72 (long int) p, mpfr_print_rnd_mode (rnd));
73 printf ("a="); mpfr_dump (a);
74 printf ("b="); mpfr_dump (b);
75 printf ("c="); mpfr_dump (c);
76 printf ("d="); mpfr_dump (d);
77 printf ("Expected\n ");
78 mpfr_dump (ref);
79 printf (" with inex = %d\n", inex_ref);
80 printf ("Got\n ");
81 mpfr_dump (res);
82 printf (" with inex = %d\n", inex_res);
83 exit (1);
84 }
85
86 mpfr_clear (ab);
87 mpfr_clear (cd);
88 mpfr_clear (res);
89 mpfr_clear (ref);
90 }
91
92 static void
93 random_tests (void)
94 {
95 mpfr_prec_t p;
96 int r;
97 mpfr_t a, b, c, d;
98
99 for (p = MPFR_PREC_MIN; p <= 4096; p++)
100 {
101 mpfr_inits2 (p, a, b, c, d, (mpfr_ptr) 0);
102 mpfr_urandomb (a, RANDS);
103 mpfr_urandomb (b, RANDS);
104 mpfr_urandomb (c, RANDS);
105 mpfr_urandomb (d, RANDS);
106 RND_LOOP_NO_RNDF (r)
107 random_test (a, b, c, d, (mpfr_rnd_t) r);
108 mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
109 }
110 }
111
112 static void
113 zero_tests (void)
114 {
115 mpfr_exp_t emin, emax;
116 mpfr_t a, b, c, d, res;
117 int i, r;
118
119 emin = mpfr_get_emin ();
120 emax = mpfr_get_emax ();
121 set_emin (MPFR_EMIN_MIN);
122 set_emax (MPFR_EMAX_MAX);
123
124 mpfr_inits2 (GMP_NUMB_BITS, a, b, c, d, (mpfr_ptr) 0);
125 for (i = 0; i <= 4; i++)
126 {
127 switch (i)
128 {
129 case 0: case 1:
130 mpfr_set_ui (a, i, MPFR_RNDN);
131 break;
132 case 2:
133 mpfr_setmax (a, MPFR_EMAX_MAX);
134 break;
135 case 3:
136 mpfr_setmin (a, MPFR_EMIN_MIN);
137 break;
138 case 4:
139 mpfr_setmin (a, MPFR_EMIN_MIN+1);
140 break;
141 }
142 RND_LOOP (r)
143 {
144 int j, inex;
145 mpfr_flags_t flags;
146
147 mpfr_set (b, a, MPFR_RNDN);
148 mpfr_set (c, a, MPFR_RNDN);
149 mpfr_neg (d, a, MPFR_RNDN);
150 /* We also want to test cases where the precision of the
151 result is twice the precision of each input, in case
152 add1sp.c and/or sub1sp.c could be involved. */
153 for (j = 1; j <= 2; j++)
154 {
155 mpfr_init2 (res, GMP_NUMB_BITS * j);
156 mpfr_clear_flags ();
157 inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r);
158 flags = __gmpfr_flags;
159 if (flags != 0 || inex != 0 || ! mpfr_zero_p (res) ||
160 (r == MPFR_RNDD ? MPFR_IS_POS (res) : MPFR_IS_NEG (res)))
161 {
162 printf ("Error in zero_tests on i = %d, %s\n",
163 i, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
164 printf ("Expected %c0, inex = 0\n",
165 r == MPFR_RNDD ? '-' : '+');
166 printf ("Got ");
167 if (MPFR_IS_POS (res))
168 printf ("+");
169 mpfr_out_str (stdout, 16, 0, res, MPFR_RNDN);
170 printf (", inex = %d\n", inex);
171 printf ("Expected flags:");
172 flags_out (0);
173 printf ("Got flags: ");
174 flags_out (flags);
175 exit (1);
176 }
177 mpfr_clear (res);
178 } /* j */
179 } /* r */
180 } /* i */
181 mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
182
183 set_emin (emin);
184 set_emax (emax);
185 }
186
187 static void
188 max_tests (void)
189 {
190 mpfr_exp_t emin, emax;
191 mpfr_t x, y1, y2;
192 int r;
193 int i, inex1, inex2;
194 mpfr_flags_t flags1, flags2;
195
196 emin = mpfr_get_emin ();
197 emax = mpfr_get_emax ();
198 set_emin (MPFR_EMIN_MIN);
199 set_emax (MPFR_EMAX_MAX);
200
201 mpfr_init2 (x, GMP_NUMB_BITS);
202 mpfr_setmax (x, MPFR_EMAX_MAX);
203 flags1 = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
204 RND_LOOP (r)
205 {
206 /* We also want to test cases where the precision of the
207 result is twice the precision of each input, in case
208 add1sp.c and/or sub1sp.c could be involved. */
209 for (i = 1; i <= 2; i++)
210 {
211 mpfr_inits2 (GMP_NUMB_BITS * i, y1, y2, (mpfr_ptr) 0);
212 inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r);
213 mpfr_clear_flags ();
214 inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r);
215 flags2 = __gmpfr_flags;
216 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
217 mpfr_equal_p (y1, y2)))
218 {
219 printf ("Error in max_tests for %s\n",
220 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
221 printf ("Expected ");
222 mpfr_dump (y1);
223 printf (" with inex = %d, flags =", inex1);
224 flags_out (flags1);
225 printf ("Got ");
226 mpfr_dump (y2);
227 printf (" with inex = %d, flags =", inex2);
228 flags_out (flags2);
229 exit (1);
230 }
231 mpfr_clears (y1, y2, (mpfr_ptr) 0);
232 } /* i */
233 } /* r */
234 mpfr_clear (x);
235
236 set_emin (emin);
237 set_emax (emax);
238 }
239
240 /* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't.
241 * a^2 + cd where a^2 overflows and cd doesn't affect the overflow
242 * (and cd may even underflow).
243 */
244 static void
245 overflow_tests (void)
246 {
247 mpfr_exp_t old_emax;
248 int i;
249
250 old_emax = mpfr_get_emax ();
251
252 for (i = 0; i < 40; i++)
253 {
254 mpfr_exp_t emax, exp_a;
255 mpfr_t a, k, c, d, z1, z2;
256 mpfr_rnd_t rnd;
257 mpfr_prec_t prec_a, prec_z;
258 int add = i & 1, inex, inex1, inex2;
259 mpfr_flags_t flags1, flags2;
260
261 /* In most cases, we do the test with the maximum exponent. */
262 emax = i % 8 != 0 ? MPFR_EMAX_MAX : 64 + RAND_BOOL ();
263 set_emax (emax);
264 exp_a = emax/2 + 32;
265
266 rnd = RND_RAND_NO_RNDF ();
267 prec_a = 64 + randlimb () % 100;
268 prec_z = MPFR_PREC_MIN + randlimb () % 160;
269
270 mpfr_init2 (a, prec_a);
271 mpfr_urandom (a, RANDS, MPFR_RNDU);
272 mpfr_set_exp (a, exp_a);
273
274 mpfr_init2 (k, prec_a - 32);
275 mpfr_urandom (k, RANDS, MPFR_RNDU);
276 mpfr_set_exp (k, exp_a - 32);
277
278 mpfr_init2 (c, prec_a + 1);
279 inex = mpfr_add (c, a, k, MPFR_RNDN);
280 MPFR_ASSERTN (inex == 0);
281
282 mpfr_init2 (d, prec_a);
283 inex = mpfr_sub (d, a, k, MPFR_RNDN);
284 MPFR_ASSERTN (inex == 0);
285 if (add)
286 mpfr_neg (d, d, MPFR_RNDN);
287
288 mpfr_init2 (z1, prec_z);
289 mpfr_clear_flags ();
290 inex1 = mpfr_sqr (z1, k, rnd);
291 flags1 = __gmpfr_flags;
292
293 mpfr_init2 (z2, prec_z);
294 mpfr_clear_flags ();
295 inex2 = add ?
296 mpfr_fmma (z2, a, a, c, d, rnd) :
297 mpfr_fmms (z2, a, a, c, d, rnd);
298 flags2 = __gmpfr_flags;
299
300 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
301 mpfr_equal_p (z1, z2)))
302 {
303 printf ("Error 1 in overflow_tests for %s\n",
304 mpfr_print_rnd_mode (rnd));
305 printf ("Expected ");
306 mpfr_dump (z1);
307 printf (" with inex = %d, flags =", inex1);
308 flags_out (flags1);
309 printf ("Got ");
310 mpfr_dump (z2);
311 printf (" with inex = %d, flags =", inex2);
312 flags_out (flags2);
313 exit (1);
314 }
315
316 /* c and d such that a^2 +/- cd ~= a^2 (overflow) */
317 mpfr_urandom (c, RANDS, MPFR_RNDU);
318 mpfr_set_exp (c, RAND_BOOL () ? exp_a - 2 : __gmpfr_emin);
319 if (RAND_BOOL ())
320 mpfr_neg (c, c, MPFR_RNDN);
321 mpfr_urandom (d, RANDS, MPFR_RNDU);
322 mpfr_set_exp (d, RAND_BOOL () ? exp_a - 2 : __gmpfr_emin);
323 if (RAND_BOOL ())
324 mpfr_neg (d, d, MPFR_RNDN);
325
326 mpfr_clear_flags ();
327 inex1 = mpfr_sqr (z1, a, rnd);
328 flags1 = __gmpfr_flags;
329 MPFR_ASSERTN (flags1 == (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT));
330
331 mpfr_clear_flags ();
332 inex2 = add ?
333 mpfr_fmma (z2, a, a, c, d, rnd) :
334 mpfr_fmms (z2, a, a, c, d, rnd);
335 flags2 = __gmpfr_flags;
336
337 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
338 mpfr_equal_p (z1, z2)))
339 {
340 printf ("Error 2 in overflow_tests for %s\n",
341 mpfr_print_rnd_mode (rnd));
342 printf ("Expected ");
343 mpfr_dump (z1);
344 printf (" with inex = %d, flags =", inex1);
345 flags_out (flags1);
346 printf ("Got ");
347 mpfr_dump (z2);
348 printf (" with inex = %d, flags =", inex2);
349 flags_out (flags2);
350 printf ("a="); mpfr_dump (a);
351 printf ("c="); mpfr_dump (c);
352 printf ("d="); mpfr_dump (d);
353 exit (1);
354 }
355
356 mpfr_clears (a, k, c, d, z1, z2, (mpfr_ptr) 0);
357 }
358
359 set_emax (old_emax);
360 }
361
362 /* (1/2)x + (1/2)x = x tested near underflow */
363 static void
364 half_plus_half (void)
365 {
366 mpfr_exp_t emin;
367 mpfr_t h, x1, x2, y;
368 int neg, r, i, inex;
369 mpfr_flags_t flags;
370
371 emin = mpfr_get_emin ();
372 set_emin (MPFR_EMIN_MIN);
373 mpfr_inits2 (4, h, x1, (mpfr_ptr) 0);
374 mpfr_init2 (x2, GMP_NUMB_BITS);
375 mpfr_set_ui_2exp (h, 1, -1, MPFR_RNDN);
376
377 for (mpfr_setmin (x1, __gmpfr_emin);
378 MPFR_GET_EXP (x1) < __gmpfr_emin + 2;
379 mpfr_nextabove (x1))
380 {
381 inex = mpfr_set (x2, x1, MPFR_RNDN);
382 MPFR_ASSERTN (inex == 0);
383 for (neg = 0; neg < 2; neg++)
384 {
385 RND_LOOP (r)
386 {
387 /* We also want to test cases where the precision of the
388 result is twice the precision of each input, in case
389 add1sp.c and/or sub1sp.c could be involved. */
390 for (i = 1; i <= 2; i++)
391 {
392 mpfr_init2 (y, GMP_NUMB_BITS * i);
393 mpfr_clear_flags ();
394 inex = mpfr_fmma (y, h, x2, h, x2, (mpfr_rnd_t) r);
395 flags = __gmpfr_flags;
396 if (! (flags == 0 && inex == 0 && mpfr_equal_p (y, x2)))
397 {
398 printf ("Error in half_plus_half for %s\n",
399 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
400 printf ("Expected ");
401 mpfr_dump (x2);
402 printf (" with inex = 0, flags =");
403 flags_out (0);
404 printf ("Got ");
405 mpfr_dump (y);
406 printf (" with inex = %d, flags =", inex);
407 flags_out (flags);
408 exit (1);
409 }
410 mpfr_clear (y);
411 }
412 }
413 mpfr_neg (x2, x2, MPFR_RNDN);
414 }
415 }
416
417 mpfr_clears (h, x1, x2, (mpfr_ptr) 0);
418 set_emin (emin);
419 }
420
421 /* check that result has exponent in [emin, emax]
422 (see https://sympa.inria.fr/sympa/arc/mpfr/2017-04/msg00016.html)
423 Overflow detection in sub1.c was incorrect (only for UBF cases);
424 fixed in r11414. */
425 static void
426 bug20170405 (void)
427 {
428 mpfr_t x, y, z;
429
430 mpfr_inits2 (866, x, y, z, (mpfr_ptr) 0);
431
432 mpfr_set_str_binary (x, "-0.10010101110110001111000110010100011001111011110001010100010000111110001010111110100001000000011010001000010000101110000000001100001011000110010000100111001100000101110000000001001101101101010110000110100010010111011001101101010011111000101100000010001100010000011100000000011110100010111011101011000101101011110110001010011001101110011101100001111000011000000011000010101010000101001001010000111101100001000001011110011000110010001100001101101001001010000111100101000010111001001101010011001110110001000010101001100000101010110000100100100010101011111001001100010001010110011000000001011110011000110001000100101000111010010111111110010111001101110101010010101101010100111001011100101101010111010011001000001010011001010001101000111011010010100110011001111111000011101111001010111001001011011011110101101001100011010001010110011100001101100100001001100111010100010100E768635654");
433 mpfr_set_str_binary (y, "-0.11010001010111110010110101010011000010010011010011011101100100110000110101100110111010001001110101110000011101100010110100100011001101111010100011111001011100111101110101101001000101011110101101101011010100110010111111010011011100101111110011001001010101011101111100011101100001010010011000110010110101001110010001101111111001100100000101010100110011101101101010011001000110100001001100000010110010101111000110110000111011000110001000100100100101111110001111100101011100100100110111010000010110110001110010001001101000000110111000101000110101111110000110001110100010101111010110001111010111111111010011001001100110011000110010110011000101110001010001101000100010000110011101010010010011110100000111100000101100110001111010000100011111000001101111110100000011011110010100010010011111111000010110000000011010011001100110001110111111010111110000111110010110011001000010E768635576");
434 /* since emax = 1073741821, x^2-y^2 should overflow */
435 mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
436 MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);
437
438 /* same test when z has a different precision */
439 mpfr_set_prec (z, 867);
440 mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
441 MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);
442
443 mpfr_set_prec (x, 564);
444 mpfr_set_prec (y, 564);
445 mpfr_set_prec (z, 2256);
446 mpfr_set_str_binary (x, "1e-536870913");
447 mpfr_set_str_binary (y, "-1e-536870913");
448 mpfr_fmma (z, x, y, x, y, MPFR_RNDN);
449 /* we should get -0 as result */
450 MPFR_ASSERTN(mpfr_zero_p (z) && mpfr_signbit (z));
451
452 mpfr_set_prec (x, 564);
453 mpfr_set_prec (y, 564);
454 mpfr_set_prec (z, 2256);
455 mpfr_set_str_binary (x, "1.10010000111100110011001101111111101000111001011000110100110010000101000100010001000000111100010000101001011011111001111000110101111100101111001100001100011101100100011110000000011000010110111100111000100101010001011111010111011001110010001011101111001011001110110000010000011100010001010001011100100110111110101001001111001011101111110011101110101010110100010010111011111100010101111100011110111001011111101110101101101110100101111010000101011110100000000110111101000001100001000100010110100111010011011010110011100111010000101110010101111001011100110101100001100e-737194993");
456 mpfr_set_str_binary (y, "-1.00101000100001001101011110100010110011101010011011010111100110101011111100000100101000111010111101100100110010001110011011100100110110000001011001000111101111101111110101100110111000000011000001101001010100010010001110001000011010000100111001001100101111111100010101110101001101101101111010100011011110001000010000010100011000011000010110101100000111111110111001100100100001101011111011100101110111000100101010110100010011101010110010100110100111000000100111101101101000000011110000100110100100011000010011110010001010000110100011111101101101110001110001101101010e-737194903");
457 mpfr_fmma (z, x, y, x, y, MPFR_RNDN);
458 /* we should get -0 as result */
459 MPFR_ASSERTN(mpfr_zero_p (z) && mpfr_signbit (z));
460
461 mpfr_set_prec (x, 2);
462 mpfr_set_prec (y, 2);
463 mpfr_set_prec (z, 2);
464 /* (a+i*b)*(c+i*d) with:
465 a=0.10E1
466 b=0.10E-536870912
467 c=0.10E-536870912
468 d=0.10E1 */
469 mpfr_set_str_binary (x, "0.10E1"); /* x = a = d */
470 mpfr_set_str_binary (y, "0.10E-536870912"); /* y = b = c */
471 /* real part is a*c-b*d = x*y-y*x */
472 mpfr_fmms (z, x, y, y, x, MPFR_RNDN);
473 MPFR_ASSERTN(mpfr_zero_p (z) && !mpfr_signbit (z));
474 /* imaginary part is a*d+b*c = x*x+y*y */
475 mpfr_fmma (z, x, x, y, y, MPFR_RNDN);
476 MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0);
477 mpfr_fmma (z, y, y, x, x, MPFR_RNDN);
478 MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0);
479
480 mpfr_clears (x, y, z, (mpfr_ptr) 0);
481 }
482
483 static void
484 underflow_tests (void)
485 {
486 mpfr_t x, y, z;
487 mpfr_prec_t p;
488 mpfr_exp_t emin;
489
490 emin = mpfr_get_emin ();
491 set_emin (-17);
492 for (p = MPFR_PREC_MIN; p <= 1024; p++)
493 {
494 mpfr_inits2 (p, x, y, (mpfr_ptr) 0);
495 mpfr_init2 (z, p + 1);
496 mpfr_set_str_binary (x, "1e-10");
497 mpfr_set_str_binary (y, "1e-11");
498 mpfr_clear_underflow ();
499 mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
500 /* the exact result is 2^-20-2^-22, and should underflow */
501 MPFR_ASSERTN(mpfr_underflow_p ());
502 MPFR_ASSERTN(mpfr_zero_p (z));
503 MPFR_ASSERTN(mpfr_signbit (z) == 0);
504 mpfr_clears (x, y, z, (mpfr_ptr) 0);
505 }
506 set_emin (emin);
507 }
508
509 static void
510 bug20180604 (void)
511 {
512 mpfr_t x, y, yneg, z;
513 mpfr_exp_t emin;
514 int ret;
515
516 emin = mpfr_get_emin ();
517 set_emin (-1073741821);
518 mpfr_inits2 (564, x, y, yneg, (mpfr_ptr) 0);
519 mpfr_init2 (z, 2256);
520 mpfr_set_str_binary (x, "1.10010000111100110011001101111111101000111001011000110100110010000101000100010001000000111100010000101001011011111001111000110101111100101111001100001100011101100100011110000000011000010110111100111000100101010001011111010111011001110010001011101111001011001110110000010000011100010001010001011100100110111110101001001111001011101111110011101110101010110100010010111011111100010101111100011110111001011111101110101101101110100101111010000101011110100000000110111101000001100001000100010110100111010011011010110011100111010000101110010101111001011100110101100001100E-737194993");
521 mpfr_set_str_binary (y, "-1.00101000100001001101011110100010110011101010011011010111100110101011111100000100101000111010111101100100110010001110011011100100110110000001011001000111101111101111110101100110111000000011000001101001010100010010001110001000011010000100111001001100101111111100010101110101001101101101111010100011011110001000010000010100011000011000010110101100000111111110111001100100100001101011111011100101110111000100101010110100010011101010110010100110100111000000100111101101101000000011110000100110100100011000010011110010001010000110100011111101101101110001110001101101010E-737194903");
522
523 mpfr_clear_underflow ();
524 ret = mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
525 MPFR_ASSERTN(mpfr_underflow_p ());
526 MPFR_ASSERTN(mpfr_zero_p (z));
527 MPFR_ASSERTN(mpfr_signbit (z) == 1);
528 MPFR_ASSERTN(ret > 0);
529
530 mpfr_clear_underflow ();
531 ret = mpfr_fmms (z, y, y, x, x, MPFR_RNDN);
532 MPFR_ASSERTN(mpfr_underflow_p ());
533 MPFR_ASSERTN(mpfr_zero_p (z));
534 MPFR_ASSERTN(mpfr_signbit (z) == 0);
535 MPFR_ASSERTN(ret < 0);
536
537 mpfr_neg (yneg, y, MPFR_RNDN);
538 mpfr_clear_underflow ();
539 ret = mpfr_fmms (z, x, x, y, yneg, MPFR_RNDN);
540 MPFR_ASSERTN(mpfr_underflow_p ());
541 MPFR_ASSERTN(mpfr_zero_p (z));
542 MPFR_ASSERTN(mpfr_signbit (z) == 0);
543 MPFR_ASSERTN(ret < 0);
544
545 mpfr_clear_underflow ();
546 ret = mpfr_fmms (z, y, yneg, x, x, MPFR_RNDN);
547 MPFR_ASSERTN(mpfr_underflow_p ());
548 MPFR_ASSERTN(mpfr_zero_p (z));
549 MPFR_ASSERTN(mpfr_signbit (z) == 1);
550 MPFR_ASSERTN(ret > 0);
551
552 mpfr_clears (x, y, yneg, z, (mpfr_ptr) 0);
553 set_emin (emin);
554 }
555
556 /* this bug was discovered from mpc_mul */
557 static void
558 bug20200206 (void)
559 {
560 mpfr_exp_t emin = mpfr_get_emin ();
561 mpfr_t xre, xim, yre, yim, zre;
562
563 set_emin (-1073);
564 mpfr_inits2 (53, xre, xim, yre, yim, zre, (mpfr_ptr) 0);
565 mpfr_set_str (xre, "-0x8.294611b331c8p-904", 16, MPFR_RNDN);
566 mpfr_set_str (xim, "-0x1.859278c2992fap-676", 16, MPFR_RNDN);
567 mpfr_set_str (yre, "0x9.ac54802a95f8p-820", 16, MPFR_RNDN);
568 mpfr_set_str (yim, "0x3.17e59e7612aap-412", 16, MPFR_RNDN);
569 mpfr_fmms (zre, xre, yre, xim, yim, MPFR_RNDN);
570 if (mpfr_regular_p (zre) && mpfr_get_exp (zre) < -1073)
571 {
572 printf ("Error, mpfr_fmms returns an out-of-range exponent:\n");
573 mpfr_dump (zre);
574 exit (1);
575 }
576 mpfr_clears (xre, xim, yre, yim, zre, (mpfr_ptr) 0);
577 set_emin (emin);
578 }
579
580 /* check for integer overflow (see bug fixed in r13798) */
581 static void
582 extreme_underflow (void)
583 {
584 mpfr_t x, y, z;
585 mpfr_exp_t emin = mpfr_get_emin ();
586
587 set_emin (MPFR_EMIN_MIN);
588 mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0);
589 mpfr_set_ui_2exp (x, 1, MPFR_EMIN_MIN - 1, MPFR_RNDN);
590 mpfr_set (y, x, MPFR_RNDN);
591 mpfr_nextabove (x);
592 mpfr_clear_flags ();
593 mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
594 MPFR_ASSERTN (__gmpfr_flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT));
595 MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z));
596 mpfr_clears (x, y, z, (mpfr_ptr) 0);
597 set_emin (emin);
598 }
599
600 /* Test double-rounding cases in mpfr_set_1_2, which is called when
601 all the precisions are the same. With set.c before r13347, this
602 triggers errors for neg=0. */
603 static void
604 double_rounding (void)
605 {
606 int p;
607
608 for (p = 4; p < 4 * GMP_NUMB_BITS; p++)
609 {
610 mpfr_t a, b, c, d, z1, z2;
611 int neg;
612
613 mpfr_inits2 (p, a, b, c, d, z1, z2, (mpfr_ptr) 0);
614 /* Choose a, b, c, d such that a*b = 2^p and c*d = 1 + epsilon */
615 mpfr_set_ui (a, 2, MPFR_RNDN);
616 mpfr_set_ui_2exp (b, 1, p - 1, MPFR_RNDN);
617 mpfr_set_ui (c, 1, MPFR_RNDN);
618 mpfr_nextabove (c);
619 /* c = 1 + ulp(1) = 1 + 2 * ulp(1/2) */
620
621 for (neg = 0; neg <= 1; neg++)
622 {
623 int inex1, inex2;
624
625 /* Set d = 1 - (1 + neg) * ulp(1/2), thus
626 * c*d = 1 + (1 - neg) * ulp(1/2) - 2 * (1 + neg) * ulp(1/2)^2,
627 * so that a*b + c*d rounds to 2^p + 1 and epsilon has the
628 * same sign as (-1)^neg.
629 */
630 mpfr_set_ui (d, 1, MPFR_RNDN);
631 mpfr_nextbelow (d);
632 mpfr_set_ui_2exp (z1, 1, p, MPFR_RNDN);
633 if (neg)
634 {
635 mpfr_nextbelow (d);
636 inex1 = -1;
637 }
638 else
639 {
640 mpfr_nextabove (z1);
641 inex1 = 1;
642 }
643
644 inex2 = mpfr_fmma (z2, a, b, c, d, MPFR_RNDN);
645 if (!(mpfr_equal_p (z1, z2) && SAME_SIGN (inex1, inex2)))
646 {
647 printf ("Error in double_rounding for precision %d, neg=%d\n",
648 p, neg);
649 printf ("Expected ");
650 mpfr_dump (z1);
651 printf ("with inex = %d\n", inex1);
652 printf ("Got ");
653 mpfr_dump (z2);
654 printf ("with inex = %d\n", inex2);
655 exit (1);
656 }
657 }
658
659 mpfr_clears (a, b, c, d, z1, z2, (mpfr_ptr) 0);
660 }
661 }
662
663 int
664 main (int argc, char *argv[])
665 {
666 tests_start_mpfr ();
667
668 bug20200206 ();
669 bug20180604 ();
670 underflow_tests ();
671 random_tests ();
672 zero_tests ();
673 max_tests ();
674 overflow_tests ();
675 half_plus_half ();
676 bug20170405 ();
677 double_rounding ();
678 extreme_underflow ();
679
680 tests_end_mpfr ();
681 return 0;
682 }