1 /* Test file for mpfr_fma.
2
3 Copyright 2001-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 /* When a * b is exact, the FMA is equivalent to the separate operations. */
26 static void
test_exact(void)27 test_exact (void)
28 {
29 const char *val[] =
30 { "@NaN@", "-@Inf@", "-2", "-1", "-0", "0", "1", "2", "@Inf@" };
31 int sv = numberof (val);
32 int i, j, k;
33 int rnd;
34 mpfr_t a, b, c, r1, r2;
35
36 mpfr_inits2 (8, a, b, c, r1, r2, (mpfr_ptr) 0);
37
38 for (i = 0; i < sv; i++)
39 for (j = 0; j < sv; j++)
40 for (k = 0; k < sv; k++)
41 RND_LOOP (rnd)
42 {
43 if (mpfr_set_str (a, val[i], 10, MPFR_RNDN) ||
44 mpfr_set_str (b, val[j], 10, MPFR_RNDN) ||
45 mpfr_set_str (c, val[k], 10, MPFR_RNDN) ||
46 mpfr_mul (r1, a, b, (mpfr_rnd_t) rnd) ||
47 mpfr_add (r1, r1, c, (mpfr_rnd_t) rnd))
48 {
49 if (rnd == MPFR_RNDF)
50 continue;
51 printf ("test_exact internal error for (%d,%d,%d,%d,%s)\n",
52 i, j, k, rnd, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
53 exit (1);
54 }
55 if (mpfr_fma (r2, a, b, c, (mpfr_rnd_t) rnd))
56 {
57 printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be exact\n",
58 i, j, k, rnd);
59 exit (1);
60 }
61 if (MPFR_IS_NAN (r1))
62 {
63 if (MPFR_IS_NAN (r2))
64 continue;
65 printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be NaN\n",
66 i, j, k, rnd);
67 exit (1);
68 }
69 if (! mpfr_equal_p (r1, r2) || MPFR_SIGN (r1) != MPFR_SIGN (r2))
70 {
71 printf ("test_exact(%d,%d,%d,%d):\nexpected ", i, j, k, rnd);
72 mpfr_out_str (stdout, 10, 0, r1, MPFR_RNDN);
73 printf ("\n got ");
74 mpfr_out_str (stdout, 10, 0, r2, MPFR_RNDN);
75 printf ("\n");
76 exit (1);
77 }
78 }
79
80 mpfr_clears (a, b, c, r1, r2, (mpfr_ptr) 0);
81 }
82
83 static void
test_overflow1(void)84 test_overflow1 (void)
85 {
86 mpfr_t x, y, z, r;
87 int inex;
88
89 mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0);
90 MPFR_SET_POS (x);
91 mpfr_setmax (x, mpfr_get_emax ()); /* x = 2^emax - ulp */
92 mpfr_set_ui (y, 2, MPFR_RNDN); /* y = 2 */
93 mpfr_neg (z, x, MPFR_RNDN); /* z = -x = -(2^emax - ulp) */
94 mpfr_clear_flags ();
95 /* The intermediate multiplication x * y overflows, but x * y + z = x
96 is representable. */
97 inex = mpfr_fma (r, x, y, z, MPFR_RNDN);
98 if (inex || ! mpfr_equal_p (r, x))
99 {
100 printf ("Error in test_overflow1\nexpected ");
101 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
102 printf (" with inex = 0\n got ");
103 mpfr_out_str (stdout, 2, 0, r, MPFR_RNDN);
104 printf (" with inex = %d\n", inex);
105 exit (1);
106 }
107 if (mpfr_overflow_p ())
108 {
109 printf ("Error in test_overflow1: overflow flag set\n");
110 exit (1);
111 }
112 mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
113 }
114
115 static void
test_overflow2(void)116 test_overflow2 (void)
117 {
118 mpfr_t x, y, z, r;
119 int i, inex, rnd, err = 0;
120
121 mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0);
122
123 MPFR_SET_POS (x);
124 mpfr_setmin (x, mpfr_get_emax ()); /* x = 0.1@emax */
125 mpfr_set_si (y, -2, MPFR_RNDN); /* y = -2 */
126 /* The intermediate multiplication x * y will overflow. */
127
128 for (i = -9; i <= 9; i++)
129 RND_LOOP_NO_RNDF (rnd)
130 {
131 int inf, overflow;
132
133 inf = rnd == MPFR_RNDN || rnd == MPFR_RNDD || rnd == MPFR_RNDA;
134 overflow = inf || i <= 0;
135
136 inex = mpfr_set_si_2exp (z, i, mpfr_get_emin (), MPFR_RNDN);
137 MPFR_ASSERTN (inex == 0);
138
139 mpfr_clear_flags ();
140 /* One has: x * y = -1@emax exactly (but not representable). */
141 inex = mpfr_fma (r, x, y, z, (mpfr_rnd_t) rnd);
142 if (overflow ^ (mpfr_overflow_p () != 0))
143 {
144 printf ("Error in test_overflow2 (i = %d, %s): wrong overflow"
145 " flag (should be %d)\n", i,
146 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), overflow);
147 err = 1;
148 }
149 if (mpfr_nanflag_p ())
150 {
151 printf ("Error in test_overflow2 (i = %d, %s): NaN flag should"
152 " not be set\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
153 err = 1;
154 }
155 if (mpfr_nan_p (r))
156 {
157 printf ("Error in test_overflow2 (i = %d, %s): got NaN\n",
158 i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
159 err = 1;
160 }
161 else if (MPFR_IS_POS (r))
162 {
163 printf ("Error in test_overflow2 (i = %d, %s): wrong sign "
164 "(+ instead of -)\n", i,
165 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
166 err = 1;
167 }
168 else if (inf && ! mpfr_inf_p (r))
169 {
170 printf ("Error in test_overflow2 (i = %d, %s): expected -Inf,"
171 " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
172 mpfr_dump (r);
173 err = 1;
174 }
175 else if (!inf && (mpfr_inf_p (r) ||
176 (mpfr_nextbelow (r), ! mpfr_inf_p (r))))
177 {
178 printf ("Error in test_overflow2 (i = %d, %s): expected -MAX,"
179 " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
180 mpfr_dump (r);
181 err = 1;
182 }
183 if (inf ? inex >= 0 : inex <= 0)
184 {
185 printf ("Error in test_overflow2 (i = %d, %s): wrong inexact"
186 " flag (got %d)\n", i,
187 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex);
188 err = 1;
189 }
190
191 }
192
193 if (err)
194 exit (1);
195 mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
196 }
197
198 static void
test_overflow3(void)199 test_overflow3 (void)
200 {
201 mpfr_t x, y, z, r;
202 int inex;
203 mpfr_prec_t p = 8;
204 mpfr_flags_t ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT, flags;
205 int i, j, k;
206 unsigned int neg;
207
208 mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0);
209 for (i = 0; i < 2; i++)
210 {
211 mpfr_init2 (r, 2 * p + i);
212 mpfr_set_ui_2exp (x, 1, mpfr_get_emax () - 1, MPFR_RNDN);
213 mpfr_set_ui (y, 2, MPFR_RNDN); /* y = 2 */
214 for (j = 1; j <= 2; j++)
215 for (k = 0; k <= 1; k++)
216 {
217 mpfr_set_si_2exp (z, -1, mpfr_get_emax () - mpfr_get_prec (r) - j,
218 MPFR_RNDN);
219 if (k)
220 mpfr_nextabove (z);
221 for (neg = 0; neg <= 3; neg++)
222 {
223 mpfr_clear_flags ();
224 /* (The following applies for neg = 0 or 2, all the signs
225 need to be reversed for neg = 1 or 3.)
226 We have x*y = 2^emax and
227 z = - 2^(emax-2p-i-j) * (1-k*2^(-p)), thus
228 x*y+z = 2^emax - 2^(emax-2p-i-j) + k*2^(emax-3p-i-j)
229 should overflow. Indeed it is >= the midpoint of
230 2^emax - 2^(emax-2p-i) and 2^emax, the midpoint
231 being obtained for j = 1 and k = 0. */
232 inex = mpfr_fma (r, x, y, z, MPFR_RNDN);
233 flags = __gmpfr_flags;
234 if (! mpfr_inf_p (r) || flags != ex_flags ||
235 ((neg & 1) == 0 ?
236 (inex <= 0 || MPFR_IS_NEG (r)) :
237 (inex >= 0 || MPFR_IS_POS (r))))
238 {
239 printf ("Error in test_overflow3 for "
240 "i=%d j=%d k=%d neg=%u\n", i, j, k, neg);
241 printf ("Expected %c@Inf@\n with inex %c 0 and flags:",
242 (neg & 1) == 0 ? '+' : '-',
243 (neg & 1) == 0 ? '>' : '<');
244 flags_out (ex_flags);
245 printf ("Got ");
246 mpfr_dump (r);
247 printf (" with inex = %d and flags:", inex);
248 flags_out (flags);
249 exit (1);
250 }
251 if (neg == 0 || neg == 2)
252 mpfr_neg (x, x, MPFR_RNDN);
253 if (neg == 1 || neg == 3)
254 mpfr_neg (y, y, MPFR_RNDN);
255 mpfr_neg (z, z, MPFR_RNDN);
256 } /* neg */
257 } /* k */
258 mpfr_clear (r);
259 } /* i */
260 mpfr_clears (x, y, z, (mpfr_ptr) 0);
261 }
262
263 static void
test_overflow4(void)264 test_overflow4 (void)
265 {
266 mpfr_t x, y, z, r1, r2;
267 mpfr_exp_t emax, e;
268 mpfr_prec_t px;
269 mpfr_flags_t flags1, flags2;
270 int inex1, inex2;
271 int ei, i, j;
272 int below;
273 unsigned int neg;
274
275 emax = mpfr_get_emax ();
276
277 mpfr_init2 (y, MPFR_PREC_MIN);
278 mpfr_set_ui (y, 2, MPFR_RNDN); /* y = 2 */
279
280 mpfr_init2 (z, 8);
281
282 for (px = 17; px < 256; px *= 2)
283 {
284 mpfr_init2 (x, px);
285 mpfr_inits2 (px - 8, r1, r2, (mpfr_ptr) 0);
286 for (ei = 0; ei <= 1; ei++)
287 {
288 e = ei ? emax : 0;
289 mpfr_set_ui_2exp (x, 1, e - 1, MPFR_RNDN);
290 mpfr_nextabove (x); /* x = 2^(e - 1) + 2^(e - px) */
291 /* x*y = 2^e + 2^(e - px + 1), which internally overflows
292 when e = emax. */
293 for (i = -4; i <= 4; i++)
294 for (j = 2; j <= 3; j++)
295 {
296 mpfr_set_si_2exp (z, -j, e - px + i, MPFR_RNDN);
297 /* If |z| <= 2^(e - px + 1), then x*y + z >= 2^e and
298 RZ(x*y + z) = 2^e with an unbounded exponent range.
299 If |z| > 2^(e - px + 1), then RZ(x*y + z) is the
300 predecessor of 2^e (since |z| < ulp(r)/2); this
301 occurs when i > 0 and when i = 0 and j > 2 */
302 mpfr_set_ui_2exp (r1, 1, e - 1, MPFR_RNDN);
303 below = i > 0 || (i == 0 && j > 2);
304 if (below)
305 mpfr_nextbelow (r1);
306 mpfr_clear_flags ();
307 inex1 = mpfr_mul_2ui (r1, r1, 1, MPFR_RNDZ);
308 if (below || e < emax)
309 {
310 inex1 = i == 0 && j == 2 ? 0 : -1;
311 flags1 = inex1 ? MPFR_FLAGS_INEXACT : 0;
312 }
313 else
314 {
315 MPFR_ASSERTN (inex1 < 0);
316 flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
317 MPFR_ASSERTN (flags1 == __gmpfr_flags);
318 }
319 for (neg = 0; neg <= 3; neg++)
320 {
321 mpfr_clear_flags ();
322 inex2 = mpfr_fma (r2, x, y, z, MPFR_RNDZ);
323 flags2 = __gmpfr_flags;
324 if (! (mpfr_equal_p (r1, r2) &&
325 SAME_SIGN (inex1, inex2) &&
326 flags1 == flags2))
327 {
328 printf ("Error in test_overflow4 for "
329 "px=%d ei=%d i=%d j=%d neg=%u\n",
330 (int) px, ei, i, j, neg);
331 printf ("Expected ");
332 mpfr_dump (r1);
333 printf ("with inex = %d and flags:", inex1);
334 flags_out (flags1);
335 printf ("Got ");
336 mpfr_dump (r2);
337 printf ("with inex = %d and flags:", inex2);
338 flags_out (flags2);
339 exit (1);
340 }
341 if (neg == 0 || neg == 2)
342 mpfr_neg (x, x, MPFR_RNDN);
343 if (neg == 1 || neg == 3)
344 mpfr_neg (y, y, MPFR_RNDN);
345 mpfr_neg (z, z, MPFR_RNDN);
346 mpfr_neg (r1, r1, MPFR_RNDN);
347 inex1 = - inex1;
348 }
349 }
350 }
351 mpfr_clears (x, r1, r2, (mpfr_ptr) 0);
352 }
353
354 mpfr_clears (y, z, (mpfr_ptr) 0);
355 }
356
357 static void
test_overflow5(void)358 test_overflow5 (void)
359 {
360 mpfr_t x, y, z, r1, r2;
361 mpfr_exp_t emax;
362 int inex1, inex2;
363 int i, rnd;
364 unsigned int neg, negr;
365
366 emax = mpfr_get_emax ();
367
368 mpfr_init2 (x, 123);
369 mpfr_init2 (y, 45);
370 mpfr_init2 (z, 67);
371 mpfr_inits2 (89, r1, r2, (mpfr_ptr) 0);
372
373 mpfr_set_ui_2exp (x, 1, emax - 1, MPFR_RNDN);
374
375 for (i = 3; i <= 17; i++)
376 {
377 mpfr_set_ui (y, i, MPFR_RNDN);
378 mpfr_set_ui_2exp (z, 1, emax - 1, MPFR_RNDN);
379 for (neg = 0; neg < 8; neg++)
380 {
381 mpfr_setsign (x, x, neg & 1, MPFR_RNDN);
382 mpfr_setsign (y, y, neg & 2, MPFR_RNDN);
383 mpfr_setsign (z, z, neg & 4, MPFR_RNDN);
384
385 /* |x*y + z| = (i +/- 1) * 2^(emax - 1) >= 2^emax (overflow)
386 and x*y + z has the same sign as x*y. */
387 negr = (neg ^ (neg >> 1)) & 1;
388
389 RND_LOOP (rnd)
390 {
391 mpfr_set_inf (r1, 1);
392 if (MPFR_IS_LIKE_RNDZ ((mpfr_rnd_t) rnd, negr))
393 {
394 mpfr_nextbelow (r1);
395 inex1 = -1;
396 }
397 else
398 inex1 = 1;
399
400 if (negr)
401 {
402 mpfr_neg (r1, r1, MPFR_RNDN);
403 inex1 = - inex1;
404 }
405
406 mpfr_clear_flags ();
407 inex2 = mpfr_fma (r2, x, y, z, (mpfr_rnd_t) rnd);
408 MPFR_ASSERTN (__gmpfr_flags ==
409 (MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW));
410
411 if (! (mpfr_equal_p (r1, r2) && SAME_SIGN (inex1, inex2)))
412 {
413 printf ("Error in test_overflow5 for i=%d neg=%u %s\n",
414 i, neg, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
415 printf ("Expected ");
416 mpfr_dump (r1);
417 printf ("with inex = %d\n", inex1);
418 printf ("Got ");
419 mpfr_dump (r2);
420 printf ("with inex = %d\n", inex2);
421 exit (1);
422 }
423 } /* rnd */
424 } /* neg */
425 } /* i */
426
427 mpfr_clears (x, y, z, r1, r2, (mpfr_ptr) 0);
428 }
429
430 static void
test_underflow1(void)431 test_underflow1 (void)
432 {
433 mpfr_t x, y, z, r;
434 int inex, signy, signz, rnd, err = 0;
435
436 mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0);
437
438 MPFR_SET_POS (x);
439 mpfr_setmin (x, mpfr_get_emin ()); /* x = 0.1@emin */
440
441 for (signy = -1; signy <= 1; signy += 2)
442 {
443 mpfr_set_si_2exp (y, signy, -1, MPFR_RNDN); /* |y| = 1/2 */
444 for (signz = -3; signz <= 3; signz += 2)
445 {
446 RND_LOOP (rnd)
447 {
448 mpfr_set_si (z, signz, MPFR_RNDN);
449 if (ABS (signz) != 1)
450 mpfr_setmax (z, mpfr_get_emax ());
451 /* |z| = 1 or 2^emax - ulp */
452 mpfr_clear_flags ();
453 inex = mpfr_fma (r, x, y, z, (mpfr_rnd_t) rnd);
454 #define STRTU1 "Error in test_underflow1 (signy = %d, signz = %d, %s)\n "
455 if (mpfr_nanflag_p ())
456 {
457 printf (STRTU1 "NaN flag is set\n", signy, signz,
458 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
459 err = 1;
460 }
461 if (signy < 0 && MPFR_IS_LIKE_RNDD(rnd, signz))
462 mpfr_nextbelow (z);
463 if (signy > 0 && MPFR_IS_LIKE_RNDU(rnd, signz))
464 mpfr_nextabove (z);
465 if ((mpfr_overflow_p () != 0) ^ (mpfr_inf_p (z) != 0))
466 {
467 printf (STRTU1 "wrong overflow flag\n", signy, signz,
468 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
469 err = 1;
470 }
471 if (mpfr_underflow_p ())
472 {
473 printf (STRTU1 "underflow flag is set\n", signy, signz,
474 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
475 err = 1;
476 }
477 if (! mpfr_equal_p (r, z))
478 {
479 printf (STRTU1 "got ", signy, signz,
480 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
481 mpfr_dump (r);
482 printf (" instead of ");
483 mpfr_dump (z);
484 err = 1;
485 }
486 if (inex >= 0 && (rnd == MPFR_RNDD ||
487 (rnd == MPFR_RNDZ && signz > 0) ||
488 (rnd == MPFR_RNDN && signy > 0)))
489 {
490 printf (STRTU1 "ternary value = %d instead of < 0\n",
491 signy, signz, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
492 inex);
493 err = 1;
494 }
495 if (inex <= 0 && (rnd == MPFR_RNDU ||
496 (rnd == MPFR_RNDZ && signz < 0) ||
497 (rnd == MPFR_RNDN && signy < 0)))
498 {
499 printf (STRTU1 "ternary value = %d instead of > 0\n",
500 signy, signz, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
501 inex);
502 err = 1;
503 }
504 }
505 }
506 }
507
508 if (err)
509 exit (1);
510 mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
511 }
512
513 static void
test_underflow2(void)514 test_underflow2 (void)
515 {
516 mpfr_t x, y, z, r1, r2;
517 int e, b, i, prec = 32, pz, inex;
518 unsigned int neg;
519
520 mpfr_init2 (x, MPFR_PREC_MIN);
521 mpfr_inits2 (prec, y, z, r1, r2, (mpfr_ptr) 0);
522
523 mpfr_set_si_2exp (x, 1, mpfr_get_emin () - 1, MPFR_RNDN);
524 /* x = 2^(emin-1) */
525
526 for (e = -1; e <= prec + 2; e++)
527 {
528 mpfr_set (z, x, MPFR_RNDN);
529 /* z = x = 2^(emin+e) */
530 for (b = 0; b <= 1; b++)
531 {
532 for (pz = prec - 4 * (b == 0); pz <= prec + 4; pz++)
533 {
534 inex = mpfr_prec_round (z, pz, MPFR_RNDZ);
535 MPFR_ASSERTN (inex == 0);
536 for (i = 15; i <= 17; i++)
537 {
538 mpfr_flags_t flags1, flags2;
539 int inex1, inex2;
540
541 mpfr_set_si_2exp (y, i, -4 - prec, MPFR_RNDN);
542
543 /* <--- r --->
544 * z = 1.000...00b with b = 0 or 1
545 * xy = 01111 (i = 15)
546 * or 10000 (i = 16)
547 * or 10001 (i = 17)
548 *
549 * x, y, and z will be modified to test the different sign
550 * combinations. In the case b = 0 (i.e. |z| is a power of
551 * 2) and x*y has a different sign from z, then y will be
552 * divided by 2, so that i = 16 corresponds to a midpoint.
553 */
554
555 for (neg = 0; neg < 8; neg++)
556 {
557 int xyneg, prev_binade;
558
559 mpfr_setsign (x, x, neg & 1, MPFR_RNDN);
560 mpfr_setsign (y, y, neg & 2, MPFR_RNDN);
561 mpfr_setsign (z, z, neg & 4, MPFR_RNDN);
562
563 xyneg = (neg ^ (neg >> 1)) & 1; /* true iff x*y < 0 */
564
565 /* If a change of binade occurs by adding x*y to z
566 exactly, then take into account the fact that the
567 midpoint has a different exponent. */
568 prev_binade = b == 0 && (xyneg ^ MPFR_IS_NEG (z));
569 if (prev_binade)
570 mpfr_div_2ui (y, y, 1, MPFR_RNDN);
571
572 mpfr_set (r1, z, MPFR_RNDN);
573 flags1 = MPFR_FLAGS_INEXACT;
574
575 if (e == -1 && i == 17 && b == 0 &&
576 (xyneg ^ (neg >> 2)) != 0)
577 {
578 /* Special underflow case. */
579 flags1 |= MPFR_FLAGS_UNDERFLOW;
580 inex1 = xyneg ? 1 : -1;
581 }
582 else if (i == 15 || (i == 16 && b == 0))
583 {
584 /* round toward z */
585 inex1 = xyneg ? 1 : -1;
586 }
587 else if (xyneg)
588 {
589 /* round away from z, with x*y < 0 */
590 mpfr_nextbelow (r1);
591 inex1 = -1;
592 }
593 else
594 {
595 /* round away from z, with x*y > 0 */
596 mpfr_nextabove (r1);
597 inex1 = 1;
598 }
599
600 mpfr_clear_flags ();
601 inex2 = mpfr_fma (r2, x, y, z, MPFR_RNDN);
602 flags2 = __gmpfr_flags;
603
604 if (! (mpfr_equal_p (r1, r2) &&
605 SAME_SIGN (inex1, inex2) &&
606 flags1 == flags2))
607 {
608 printf ("Error in test_underflow2 for "
609 "e=%d b=%d pz=%d i=%d neg=%u\n",
610 e, b, pz, i, neg);
611 printf ("Expected ");
612 mpfr_dump (r1);
613 printf ("with inex = %d and flags:", inex1);
614 flags_out (flags1);
615 printf ("Got ");
616 mpfr_dump (r2);
617 printf ("with inex = %d and flags:", inex2);
618 flags_out (flags2);
619 exit (1);
620 }
621
622 /* Restore y. */
623 if (prev_binade)
624 mpfr_mul_2ui (y, y, 1, MPFR_RNDN);
625 } /* neg */
626 } /* i */
627 } /* pz */
628
629 inex = mpfr_prec_round (z, prec, MPFR_RNDZ);
630 MPFR_ASSERTN (inex == 0);
631 MPFR_SET_POS (z);
632 mpfr_nextabove (z);
633 } /* b */
634 mpfr_mul_2ui (x, x, 1, MPFR_RNDN);
635 } /* e */
636
637 mpfr_clears (x, y, z, r1, r2, (mpfr_ptr) 0);
638 }
639
640 static void
test_underflow3(int n)641 test_underflow3 (int n)
642 {
643 mpfr_t x, y, z, t1, t2;
644 int sign, k, s, rnd, inex1, inex2;
645 mpfr_exp_t e;
646 mpfr_flags_t flags1, flags2;
647
648 mpfr_inits2 (4, x, z, t1, t2, (mpfr_ptr) 0);
649 mpfr_init2 (y, 6);
650
651 e = mpfr_get_emin () - 1;
652
653 for (sign = 1; sign >= -1; sign -= 2)
654 for (k = 1; k <= 7; k++)
655 for (s = -1; s <= 1; s++)
656 {
657 mpfr_set_si_2exp (x, sign, e, MPFR_RNDN);
658 mpfr_set_si_2exp (y, 8*k+s, -6, MPFR_RNDN);
659 mpfr_neg (z, x, MPFR_RNDN);
660 /* x = sign * 2^(emin-1)
661 y = (8 * k + s) * 2^(-6) = k / 8 + s * 2^(-6)
662 z = -x = -sign * 2^(emin-1)
663 FMA(x,y,z) = sign * ((k-8) * 2^(emin-4) + s * 2^(emin-7)) exactly.
664 Note: The purpose of the s * 2^(emin-7) term is to yield
665 double rounding when scaling for k = 4, s != 0, MPFR_RNDN. */
666
667 RND_LOOP_NO_RNDF (rnd)
668 {
669 mpfr_clear_flags ();
670 inex1 = mpfr_set_si_2exp (t1, sign * (8*k+s-64), e-6,
671 (mpfr_rnd_t) rnd);
672 flags1 = __gmpfr_flags;
673
674 mpfr_clear_flags ();
675 inex2 = mpfr_fma (t2, x, y, z, (mpfr_rnd_t) rnd);
676 flags2 = __gmpfr_flags;
677
678 if (! (mpfr_equal_p (t1, t2) &&
679 SAME_SIGN (inex1, inex2) &&
680 flags1 == flags2))
681 {
682 printf ("Error in test_underflow3, n = %d, sign = %d,"
683 " k = %d, s = %d, %s\n", n, sign, k, s,
684 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
685 printf ("Expected ");
686 mpfr_dump (t1);
687 printf (" with inex ~ %d, flags =", inex1);
688 flags_out (flags1);
689 printf ("Got ");
690 mpfr_dump (t2);
691 printf (" with inex ~ %d, flags =", inex2);
692 flags_out (flags2);
693 exit (1);
694 }
695 }
696 }
697
698 mpfr_clears (x, y, z, t1, t2, (mpfr_ptr) 0);
699 }
700
701 /* Test s = x*y + z with PREC(z) > PREC(s) + 1, x*y underflows, where
702 z + x*y and z + sign(x*y) * 2^(emin-1) do not give the same result.
703 x = 2^emin
704 y = 2^(-8)
705 z = 2^emin * (2^PREC(s) + k - 2^(-1))
706 with k = 3 for MPFR_RNDN and k = 2 for the directed rounding modes.
707 Also test the opposite versions with neg != 0.
708 */
709 static void
test_underflow4(void)710 test_underflow4 (void)
711 {
712 mpfr_t x, y, z, s1, s2;
713 mpfr_prec_t ps = 32;
714 int inex, rnd;
715
716 mpfr_inits2 (MPFR_PREC_MIN, x, y, (mpfr_ptr) 0);
717 mpfr_inits2 (ps, s1, s2, (mpfr_ptr) 0);
718 mpfr_init2 (z, ps + 2);
719
720 inex = mpfr_set_si_2exp (x, 1, mpfr_get_emin (), MPFR_RNDN);
721 MPFR_ASSERTN (inex == 0);
722 inex = mpfr_set_si_2exp (y, 1, -8, MPFR_RNDN);
723 MPFR_ASSERTN (inex == 0);
724
725 RND_LOOP_NO_RNDF (rnd)
726 {
727 mpfr_flags_t flags1, flags2;
728 int inex1, inex2;
729 unsigned int neg;
730
731 inex = mpfr_set_si_2exp (z, 1 << 1, ps, MPFR_RNDN);
732 MPFR_ASSERTN (inex == 0);
733 inex = mpfr_sub_ui (z, z, 1, MPFR_RNDN);
734 MPFR_ASSERTN (inex == 0);
735 inex = mpfr_div_2ui (z, z, 1, MPFR_RNDN);
736 MPFR_ASSERTN (inex == 0);
737 inex = mpfr_add_ui (z, z, rnd == MPFR_RNDN ? 3 : 2, MPFR_RNDN);
738 MPFR_ASSERTN (inex == 0);
739 inex = mpfr_mul (z, z, x, MPFR_RNDN);
740 MPFR_ASSERTN (inex == 0);
741
742 for (neg = 0; neg <= 3; neg++)
743 {
744 inex1 = mpfr_set (s1, z, (mpfr_rnd_t) rnd);
745 flags1 = MPFR_FLAGS_INEXACT;
746
747 mpfr_clear_flags ();
748 inex2 = mpfr_fma (s2, x, y, z, (mpfr_rnd_t) rnd);
749 flags2 = __gmpfr_flags;
750
751 if (! (mpfr_equal_p (s1, s2) &&
752 SAME_SIGN (inex1, inex2) &&
753 flags1 == flags2))
754 {
755 printf ("Error in test_underflow4 for neg=%u %s\n",
756 neg, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
757 printf ("Expected ");
758 mpfr_dump (s1);
759 printf (" with inex ~ %d, flags =", inex1);
760 flags_out (flags1);
761 printf ("Got ");
762 mpfr_dump (s2);
763 printf (" with inex ~ %d, flags =", inex2);
764 flags_out (flags2);
765 exit (1);
766 }
767
768 if (neg == 0 || neg == 2)
769 mpfr_neg (x, x, MPFR_RNDN);
770 if (neg == 1 || neg == 3)
771 mpfr_neg (y, y, MPFR_RNDN);
772 mpfr_neg (z, z, MPFR_RNDN);
773 }
774 }
775
776 mpfr_clears (x, y, z, s1, s2, (mpfr_ptr) 0);
777 }
778
779 /* Test s = x*y + z on:
780 x = +/- 2^emin
781 y = +/- 2^(-3)
782 z = +/- 2^(emin + PREC(s)) and MPFR numbers close to this value.
783 with PREC(z) from PREC(s) - 2 to PREC(s) + 8.
784 */
785 static void
test_underflow5(void)786 test_underflow5 (void)
787 {
788 mpfr_t w, x, y, z, s1, s2, t;
789 mpfr_exp_t emin;
790 mpfr_prec_t ps = 32;
791 int inex, i, j, rnd;
792 unsigned int neg;
793
794 mpfr_inits2 (MPFR_PREC_MIN, w, x, y, (mpfr_ptr) 0);
795 mpfr_inits2 (ps, s1, s2, (mpfr_ptr) 0);
796 mpfr_init2 (t, ps + 9);
797
798 emin = mpfr_get_emin ();
799
800 inex = mpfr_set_si_2exp (w, 1, emin, MPFR_RNDN); /* w = 2^emin */
801 MPFR_ASSERTN (inex == 0);
802 inex = mpfr_set_si (x, 1, MPFR_RNDN);
803 MPFR_ASSERTN (inex == 0);
804 inex = mpfr_set_si_2exp (y, 1, -3, MPFR_RNDN);
805 MPFR_ASSERTN (inex == 0);
806
807 for (i = -2; i <= 8; i++)
808 {
809 mpfr_init2 (z, ps + i);
810 inex = mpfr_set_si_2exp (z, 1, ps, MPFR_RNDN);
811 MPFR_ASSERTN (inex == 0);
812
813 for (j = 1; j <= 32; j++)
814 mpfr_nextbelow (z);
815
816 for (j = -32; j <= 32; j++)
817 {
818 for (neg = 0; neg < 8; neg++)
819 {
820 mpfr_setsign (x, x, neg & 1, MPFR_RNDN);
821 mpfr_setsign (y, y, neg & 2, MPFR_RNDN);
822 mpfr_setsign (z, z, neg & 4, MPFR_RNDN);
823
824 inex = mpfr_fma (t, x, y, z, MPFR_RNDN);
825 MPFR_ASSERTN (inex == 0);
826
827 inex = mpfr_mul (x, x, w, MPFR_RNDN);
828 MPFR_ASSERTN (inex == 0);
829 inex = mpfr_mul (z, z, w, MPFR_RNDN);
830 MPFR_ASSERTN (inex == 0);
831
832 RND_LOOP_NO_RNDF (rnd)
833 {
834 mpfr_flags_t flags1, flags2;
835 int inex1, inex2;
836
837 mpfr_clear_flags ();
838 inex1 = mpfr_mul (s1, w, t, (mpfr_rnd_t) rnd);
839 flags1 = __gmpfr_flags;
840
841 mpfr_clear_flags ();
842 inex2 = mpfr_fma (s2, x, y, z, (mpfr_rnd_t) rnd);
843 flags2 = __gmpfr_flags;
844
845 if (! (mpfr_equal_p (s1, s2) &&
846 SAME_SIGN (inex1, inex2) &&
847 flags1 == flags2))
848 {
849 printf ("Error in test_underflow5 on "
850 "i=%d j=%d neg=%u %s\n", i, j, neg,
851 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
852 printf ("Expected ");
853 mpfr_dump (s1);
854 printf (" with inex ~ %d, flags =", inex1);
855 flags_out (flags1);
856 printf ("Got ");
857 mpfr_dump (s2);
858 printf (" with inex ~ %d, flags =", inex2);
859 flags_out (flags2);
860 exit (1);
861 }
862 } /* rnd */
863
864 inex = mpfr_div (x, x, w, MPFR_RNDN);
865 MPFR_ASSERTN (inex == 0);
866 inex = mpfr_div (z, z, w, MPFR_RNDN);
867 MPFR_ASSERTN (inex == 0);
868 } /* neg */
869
870 MPFR_SET_POS (z); /* restore the value before the loop on neg */
871 mpfr_nextabove (z);
872 } /* j */
873
874 mpfr_clear (z);
875 } /* i */
876
877 mpfr_clears (w, x, y, s1, s2, t, (mpfr_ptr) 0);
878 }
879
880 static void
bug20101018(void)881 bug20101018 (void)
882 {
883 mpfr_t x, y, z, t, u;
884 int i;
885
886 mpfr_init2 (x, 64);
887 mpfr_init2 (y, 64);
888 mpfr_init2 (z, 64);
889 mpfr_init2 (t, 64);
890 mpfr_init2 (u, 64);
891
892 mpfr_set_str (x, "0xf.fffffffffffffffp-14766", 16, MPFR_RNDN);
893 mpfr_set_str (y, "-0xf.fffffffffffffffp+317", 16, MPFR_RNDN);
894 mpfr_set_str (z, "0x8.3ffffffffffe3ffp-14443", 16, MPFR_RNDN);
895 mpfr_set_str (t, "0x8.7ffffffffffc7ffp-14444", 16, MPFR_RNDN);
896 i = mpfr_fma (u, x, y, z, MPFR_RNDN);
897 if (! mpfr_equal_p (u, t))
898 {
899 printf ("Wrong result in bug20101018 (a)\n");
900 printf ("Expected ");
901 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
902 printf ("\nGot ");
903 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
904 printf ("\n");
905 exit (1);
906 }
907 if (i <= 0)
908 {
909 printf ("Wrong ternary value in bug20101018 (a)\n");
910 printf ("Expected > 0\n");
911 printf ("Got %d\n", i);
912 exit (1);
913 }
914
915 mpfr_set_str (x, "-0xf.fffffffffffffffp-11420", 16, MPFR_RNDN);
916 mpfr_set_str (y, "0xf.fffffffffffffffp+9863", 16, MPFR_RNDN);
917 mpfr_set_str (z, "0x8.fffff80ffffffffp-1551", 16, MPFR_RNDN);
918 mpfr_set_str (t, "0x9.fffff01ffffffffp-1552", 16, MPFR_RNDN);
919 i = mpfr_fma (u, x, y, z, MPFR_RNDN);
920 if (! mpfr_equal_p (u, t))
921 {
922 printf ("Wrong result in bug20101018 (b)\n");
923 printf ("Expected ");
924 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
925 printf ("\nGot ");
926 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
927 printf ("\n");
928 exit (1);
929 }
930 if (i <= 0)
931 {
932 printf ("Wrong ternary value in bug20101018 (b)\n");
933 printf ("Expected > 0\n");
934 printf ("Got %d\n", i);
935 exit (1);
936 }
937
938 mpfr_set_str (x, "0xf.fffffffffffffffp-2125", 16, MPFR_RNDN);
939 mpfr_set_str (y, "-0xf.fffffffffffffffp-6000", 16, MPFR_RNDN);
940 mpfr_set_str (z, "0x8p-8119", 16, MPFR_RNDN);
941 mpfr_set_str (t, "0x8.000000000000001p-8120", 16, MPFR_RNDN);
942 i = mpfr_fma (u, x, y, z, MPFR_RNDN);
943 if (! mpfr_equal_p (u, t))
944 {
945 printf ("Wrong result in bug20101018 (c)\n");
946 printf ("Expected ");
947 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
948 printf ("\nGot ");
949 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
950 printf ("\n");
951 exit (1);
952 }
953 if (i <= 0)
954 {
955 printf ("Wrong ternary value in bug20101018 (c)\n");
956 printf ("Expected > 0\n");
957 printf ("Got %d\n", i);
958 exit (1);
959 }
960
961 mpfr_clear (x);
962 mpfr_clear (y);
963 mpfr_clear (z);
964 mpfr_clear (t);
965 mpfr_clear (u);
966 }
967
968 /* bug found with GMP_CHECK_RANDOMIZE=1514407719 */
969 static void
bug20171219(void)970 bug20171219 (void)
971 {
972 mpfr_t x, y, z, t, u;
973 int i;
974
975 mpfr_init2 (x, 60);
976 mpfr_init2 (y, 60);
977 mpfr_init2 (z, 60);
978 mpfr_init2 (t, 68);
979 mpfr_init2 (u, 68);
980
981 mpfr_set_ui (x, 1, MPFR_RNDN);
982 mpfr_set_ui (y, 1, MPFR_RNDN);
983 mpfr_set_ui (z, 1, MPFR_RNDN);
984 mpfr_set_ui (t, 2, MPFR_RNDN);
985 i = mpfr_fma (u, x, y, z, MPFR_RNDA);
986 if (! mpfr_equal_p (u, t))
987 {
988 printf ("Wrong result in bug20171219 (a)\n");
989 printf ("Expected ");
990 mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
991 printf ("\nGot ");
992 mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
993 printf ("\n");
994 exit (1);
995 }
996 if (i != 0)
997 {
998 printf ("Wrong ternary value in bug20171219\n");
999 printf ("Expected 0\n");
1000 printf ("Got %d\n", i);
1001 exit (1);
1002 }
1003
1004 mpfr_clear (x);
1005 mpfr_clear (y);
1006 mpfr_clear (z);
1007 mpfr_clear (t);
1008 mpfr_clear (u);
1009 }
1010
1011 /* coverage test for mpfr_set_1_2, case prec < GMP_NUMB_BITS,
1012 inex > 0, rb <> 0, sb = 0 */
1013 static void
coverage(void)1014 coverage (void)
1015 {
1016 mpfr_t x, y, z, s;
1017 int inex;
1018 mpfr_exp_t emax;
1019
1020 mpfr_inits2 (GMP_NUMB_BITS - 1, x, y, z, s, (mpfr_ptr) 0);
1021
1022 /* set x to 2^(GMP_NUMB_BITS/2) + 1, y to 2^(GMP_NUMB_BITS/2) - 1 */
1023 MPFR_ASSERTN((GMP_NUMB_BITS % 2) == 0);
1024 mpfr_set_ui_2exp (x, 1, GMP_NUMB_BITS / 2, MPFR_RNDN);
1025 mpfr_sub_ui (y, x, 1, MPFR_RNDN);
1026 mpfr_add_ui (x, x, 1, MPFR_RNDN);
1027 /* we have x*y = 2^GMP_NUMB_BITS - 1, thus has exponent GMP_NUMB_BITS */
1028 /* set z to 2^(1-GMP_NUMB_BITS), with exponent 2-GMP_NUMB_BITS */
1029 mpfr_set_ui_2exp (z, 1, 1 - GMP_NUMB_BITS, MPFR_RNDN);
1030 inex = mpfr_fms (s, x, y, z, MPFR_RNDN);
1031 /* s should be 2^GMP_NUMB_BITS-2 */
1032 MPFR_ASSERTN(inex < 0);
1033 inex = mpfr_add_ui (s, s, 2, MPFR_RNDN);
1034 MPFR_ASSERTN(mpfr_cmp_ui_2exp (s, 1, GMP_NUMB_BITS) == 0);
1035
1036 /* same example, but with overflow */
1037 emax = mpfr_get_emax ();
1038 mpfr_set_emax (GMP_NUMB_BITS + 1);
1039 mpfr_set_ui_2exp (z, 1, mpfr_get_emax () - 1, MPFR_RNDZ);
1040 mpfr_clear_flags ();
1041 inex = mpfr_fma (s, x, y, z, MPFR_RNDN);
1042 /* x*y = 2^GMP_NUMB_BITS - 1, z = 2^GMP_NUMB_BITS, thus
1043 x*y+z = 2^(GMP_NUMB_BITS+1) - 1 should round to 2^(GMP_NUMB_BITS+1),
1044 thus give an overflow */
1045 MPFR_ASSERTN(inex > 0);
1046 MPFR_ASSERTN(mpfr_inf_p (s) && mpfr_sgn (s) > 0);
1047 MPFR_ASSERTN(mpfr_overflow_p ());
1048 mpfr_set_emax (emax);
1049
1050 mpfr_clear (x);
1051 mpfr_clear (y);
1052 mpfr_clear (z);
1053 mpfr_clear (s);
1054 }
1055
1056 int
main(int argc,char * argv[])1057 main (int argc, char *argv[])
1058 {
1059 mpfr_t x, y, z, s;
1060 mpfr_exp_t emin, emax;
1061 int i;
1062
1063 tests_start_mpfr ();
1064
1065 coverage ();
1066
1067 emin = mpfr_get_emin ();
1068 emax = mpfr_get_emax ();
1069
1070 bug20171219 ();
1071 bug20101018 ();
1072
1073 mpfr_init (x);
1074 mpfr_init (s);
1075 mpfr_init (y);
1076 mpfr_init (z);
1077
1078 /* check special cases */
1079 mpfr_set_prec (x, 2);
1080 mpfr_set_prec (y, 2);
1081 mpfr_set_prec (z, 2);
1082 mpfr_set_prec (s, 2);
1083 mpfr_set_str (x, "-0.75", 10, MPFR_RNDN);
1084 mpfr_set_str (y, "0.5", 10, MPFR_RNDN);
1085 mpfr_set_str (z, "0.375", 10, MPFR_RNDN);
1086 mpfr_fma (s, x, y, z, MPFR_RNDU); /* result is 0 */
1087 if (mpfr_cmp_ui (s, 0))
1088 {
1089 printf ("Error: -0.75 * 0.5 + 0.375 should be equal to 0 for prec=2\n");
1090 printf ("got instead ");
1091 mpfr_dump (s);
1092 exit (1);
1093 }
1094
1095 mpfr_set_prec (x, 27);
1096 mpfr_set_prec (y, 27);
1097 mpfr_set_prec (z, 27);
1098 mpfr_set_prec (s, 27);
1099 mpfr_set_str_binary (x, "1.11111111111111111111111111e-1");
1100 mpfr_set (y, x, MPFR_RNDN);
1101 mpfr_set_str_binary (z, "-1.00011110100011001011001001e-1");
1102 if (mpfr_fma (s, x, y, z, MPFR_RNDN) >= 0)
1103 {
1104 printf ("Wrong inexact flag for x=y=1-2^(-27)\n");
1105 exit (1);
1106 }
1107
1108 mpfr_set_nan (x);
1109 mpfr_urandomb (y, RANDS);
1110 mpfr_urandomb (z, RANDS);
1111 mpfr_fma (s, x, y, z, MPFR_RNDN);
1112 if (!mpfr_nan_p (s))
1113 {
1114 printf ("evaluation of function in x=NAN does not return NAN\n");
1115 exit (1);
1116 }
1117
1118 mpfr_set_nan (y);
1119 mpfr_urandomb (x, RANDS);
1120 mpfr_urandomb (z, RANDS);
1121 mpfr_fma (s, x, y, z, MPFR_RNDN);
1122 if (!mpfr_nan_p(s))
1123 {
1124 printf ("evaluation of function in y=NAN does not return NAN\n");
1125 exit (1);
1126 }
1127
1128 mpfr_set_nan (z);
1129 mpfr_urandomb (y, RANDS);
1130 mpfr_urandomb (x, RANDS);
1131 mpfr_fma (s, x, y, z, MPFR_RNDN);
1132 if (!mpfr_nan_p (s))
1133 {
1134 printf ("evaluation of function in z=NAN does not return NAN\n");
1135 exit (1);
1136 }
1137
1138 mpfr_set_inf (x, 1);
1139 mpfr_set_inf (y, 1);
1140 mpfr_set_inf (z, 1);
1141 mpfr_fma (s, x, y, z, MPFR_RNDN);
1142 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
1143 {
1144 printf ("Error for (+inf) * (+inf) + (+inf)\n");
1145 exit (1);
1146 }
1147
1148 mpfr_set_inf (x, -1);
1149 mpfr_set_inf (y, -1);
1150 mpfr_set_inf (z, 1);
1151 mpfr_fma (s, x, y, z, MPFR_RNDN);
1152 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
1153 {
1154 printf ("Error for (-inf) * (-inf) + (+inf)\n");
1155 exit (1);
1156 }
1157
1158 mpfr_set_inf (x, 1);
1159 mpfr_set_inf (y, -1);
1160 mpfr_set_inf (z, -1);
1161 mpfr_fma (s, x, y, z, MPFR_RNDN);
1162 if (!mpfr_inf_p (s) || mpfr_sgn (s) > 0)
1163 {
1164 printf ("Error for (+inf) * (-inf) + (-inf)\n");
1165 exit (1);
1166 }
1167
1168 mpfr_set_inf (x, -1);
1169 mpfr_set_inf (y, 1);
1170 mpfr_set_inf (z, -1);
1171 mpfr_fma (s, x, y, z, MPFR_RNDN);
1172 if (!mpfr_inf_p (s) || mpfr_sgn (s) > 0)
1173 {
1174 printf ("Error for (-inf) * (+inf) + (-inf)\n");
1175 exit (1);
1176 }
1177
1178 mpfr_set_inf (x, 1);
1179 mpfr_set_ui (y, 0, MPFR_RNDN);
1180 mpfr_urandomb (z, RANDS);
1181 mpfr_fma (s, x, y, z, MPFR_RNDN);
1182 if (!mpfr_nan_p (s))
1183 {
1184 printf ("evaluation of function in x=INF y=0 does not return NAN\n");
1185 exit (1);
1186 }
1187
1188 mpfr_set_inf (y, 1);
1189 mpfr_set_ui (x, 0, MPFR_RNDN);
1190 mpfr_urandomb (z, RANDS);
1191 mpfr_fma (s, x, y, z, MPFR_RNDN);
1192 if (!mpfr_nan_p (s))
1193 {
1194 printf ("evaluation of function in x=0 y=INF does not return NAN\n");
1195 exit (1);
1196 }
1197
1198 mpfr_set_inf (x, 1);
1199 mpfr_urandomb (y, RANDS); /* always positive */
1200 mpfr_set_inf (z, -1);
1201 mpfr_fma (s, x, y, z, MPFR_RNDN);
1202 if (!mpfr_nan_p (s))
1203 {
1204 printf ("evaluation of function in x=INF y>0 z=-INF does not return NAN\n");
1205 exit (1);
1206 }
1207
1208 mpfr_set_inf (y, 1);
1209 mpfr_urandomb (x, RANDS);
1210 mpfr_set_inf (z, -1);
1211 mpfr_fma (s, x, y, z, MPFR_RNDN);
1212 if (!mpfr_nan_p (s))
1213 {
1214 printf ("evaluation of function in x>0 y=INF z=-INF does not return NAN\n");
1215 exit (1);
1216 }
1217
1218 mpfr_set_inf (x, 1);
1219 do mpfr_urandomb (y, RANDS); while (MPFR_IS_ZERO(y));
1220 mpfr_urandomb (z, RANDS);
1221 mpfr_fma (s, x, y, z, MPFR_RNDN);
1222 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
1223 {
1224 printf ("evaluation of function in x=INF does not return INF\n");
1225 exit (1);
1226 }
1227
1228 mpfr_set_inf (y, 1);
1229 do mpfr_urandomb (x, RANDS); while (MPFR_IS_ZERO(x));
1230 mpfr_urandomb (z, RANDS);
1231 mpfr_fma (s, x, y, z, MPFR_RNDN);
1232 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
1233 {
1234 printf ("evaluation of function at y=INF does not return INF\n");
1235 exit (1);
1236 }
1237
1238 mpfr_set_inf (z, 1);
1239 mpfr_urandomb (x, RANDS);
1240 mpfr_urandomb (y, RANDS);
1241 mpfr_fma (s, x, y, z, MPFR_RNDN);
1242 if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
1243 {
1244 printf ("evaluation of function in z=INF does not return INF\n");
1245 exit (1);
1246 }
1247
1248 mpfr_set_ui (x, 0, MPFR_RNDN);
1249 mpfr_urandomb (y, RANDS);
1250 mpfr_urandomb (z, RANDS);
1251 mpfr_fma (s, x, y, z, MPFR_RNDN);
1252 if (! mpfr_equal_p (s, z))
1253 {
1254 printf ("evaluation of function in x=0 does not return z\n");
1255 exit (1);
1256 }
1257
1258 mpfr_set_ui (y, 0, MPFR_RNDN);
1259 mpfr_urandomb (x, RANDS);
1260 mpfr_urandomb (z, RANDS);
1261 mpfr_fma (s, x, y, z, MPFR_RNDN);
1262 if (! mpfr_equal_p (s, z))
1263 {
1264 printf ("evaluation of function in y=0 does not return z\n");
1265 exit (1);
1266 }
1267
1268 {
1269 mpfr_prec_t prec;
1270 mpfr_t t, slong;
1271 mpfr_rnd_t rnd;
1272 int inexact, compare;
1273 unsigned int n;
1274
1275 mpfr_prec_t p0 = 2, p1 = 200;
1276 unsigned int N = 200;
1277
1278 mpfr_init (t);
1279 mpfr_init (slong);
1280
1281 /* generic test */
1282 for (prec = p0; prec <= p1; prec++)
1283 {
1284 mpfr_set_prec (x, prec);
1285 mpfr_set_prec (y, prec);
1286 mpfr_set_prec (z, prec);
1287 mpfr_set_prec (s, prec);
1288 mpfr_set_prec (t, prec);
1289
1290 for (n = 0; n < N; n++)
1291 {
1292 mpfr_urandomb (x, RANDS);
1293 mpfr_urandomb (y, RANDS);
1294 mpfr_urandomb (z, RANDS);
1295
1296 if (randlimb () % 2)
1297 mpfr_neg (x, x, MPFR_RNDN);
1298 if (randlimb () % 2)
1299 mpfr_neg (y, y, MPFR_RNDN);
1300 if (randlimb () % 2)
1301 mpfr_neg (z, z, MPFR_RNDN);
1302
1303 rnd = RND_RAND_NO_RNDF ();
1304 mpfr_set_prec (slong, 2 * prec);
1305 if (mpfr_mul (slong, x, y, rnd))
1306 {
1307 printf ("x*y should be exact\n");
1308 exit (1);
1309 }
1310 compare = mpfr_add (t, slong, z, rnd);
1311 inexact = mpfr_fma (s, x, y, z, rnd);
1312 if (! mpfr_equal_p (s, t))
1313 {
1314 printf ("results differ for\n");
1315 printf (" x=");
1316 mpfr_dump (x);
1317 printf (" y=");
1318 mpfr_dump (y);
1319 printf (" z=");
1320 mpfr_dump (z);
1321 printf (" with prec=%u rnd_mode=%s\n", (unsigned int) prec,
1322 mpfr_print_rnd_mode (rnd));
1323 printf ("got ");
1324 mpfr_dump (s);
1325 printf ("expected ");
1326 mpfr_dump (t);
1327 printf ("approx ");
1328 mpfr_dump (slong);
1329 exit (1);
1330 }
1331 if (! SAME_SIGN (inexact, compare))
1332 {
1333 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n",
1334 mpfr_print_rnd_mode (rnd), compare, inexact);
1335 printf (" x="); mpfr_dump (x);
1336 printf (" y="); mpfr_dump (y);
1337 printf (" z="); mpfr_dump (z);
1338 printf (" s="); mpfr_dump (s);
1339 exit (1);
1340 }
1341 }
1342 }
1343 mpfr_clear (t);
1344 mpfr_clear (slong);
1345
1346 }
1347 mpfr_clear (x);
1348 mpfr_clear (y);
1349 mpfr_clear (z);
1350 mpfr_clear (s);
1351
1352 test_exact ();
1353
1354 for (i = 0; i <= 2; i++)
1355 {
1356 if (i == 0)
1357 {
1358 /* corresponds to the generic code + mpfr_check_range */
1359 set_emin (-1024);
1360 set_emax (1024);
1361 }
1362 else if (i == 1)
1363 {
1364 set_emin (MPFR_EMIN_MIN);
1365 set_emax (MPFR_EMAX_MAX);
1366 }
1367 else
1368 {
1369 MPFR_ASSERTN (i == 2);
1370 if (emin == MPFR_EMIN_MIN && emax == MPFR_EMAX_MAX)
1371 break;
1372 set_emin (emin);
1373 set_emax (emax);
1374 }
1375
1376 test_overflow1 ();
1377 test_overflow2 ();
1378 test_overflow3 ();
1379 test_overflow4 ();
1380 test_overflow5 ();
1381 test_underflow1 ();
1382 test_underflow2 ();
1383 test_underflow3 (i);
1384 test_underflow4 ();
1385 test_underflow5 ();
1386 }
1387
1388 tests_end_mpfr ();
1389 return 0;
1390 }
1391