xref: /netbsd-src/external/lgpl3/mpfr/dist/tests/tfma.c (revision a24efa7dea9f1f56c3bdb15a927d3516792ace1c)
1 /* Test file for mpfr_fma.
2 
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel 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 http://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 <stdio.h>
24 #include <stdlib.h>
25 
26 #include "mpfr-test.h"
27 
28 /* When a * b is exact, the FMA is equivalent to the separate operations. */
29 static void
30 test_exact (void)
31 {
32   const char *val[] =
33     { "@NaN@", "-@Inf@", "-2", "-1", "-0", "0", "1", "2", "@Inf@" };
34   int sv = sizeof (val) / sizeof (*val);
35   int i, j, k;
36   int rnd;
37   mpfr_t a, b, c, r1, r2;
38 
39   mpfr_inits2 (8, a, b, c, r1, r2, (mpfr_ptr) 0);
40 
41   for (i = 0; i < sv; i++)
42     for (j = 0; j < sv; j++)
43       for (k = 0; k < sv; k++)
44         RND_LOOP (rnd)
45           {
46             if (mpfr_set_str (a, val[i], 10, MPFR_RNDN) ||
47                 mpfr_set_str (b, val[j], 10, MPFR_RNDN) ||
48                 mpfr_set_str (c, val[k], 10, MPFR_RNDN) ||
49                 mpfr_mul (r1, a, b, (mpfr_rnd_t) rnd) ||
50                 mpfr_add (r1, r1, c, (mpfr_rnd_t) rnd))
51               {
52                 printf ("test_exact internal error for (%d,%d,%d,%d)\n",
53                         i, j, k, rnd);
54                 exit (1);
55               }
56             if (mpfr_fma (r2, a, b, c, (mpfr_rnd_t) rnd))
57               {
58                 printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be exact\n",
59                         i, j, k, rnd);
60                 exit (1);
61               }
62             if (MPFR_IS_NAN (r1))
63               {
64                 if (MPFR_IS_NAN (r2))
65                   continue;
66                 printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be NaN\n",
67                         i, j, k, rnd);
68                 exit (1);
69               }
70             if (mpfr_cmp (r1, r2) || MPFR_SIGN (r1) != MPFR_SIGN (r2))
71               {
72                 printf ("test_exact(%d,%d,%d,%d):\nexpected ", i, j, k, rnd);
73                 mpfr_out_str (stdout, 10, 0, r1, MPFR_RNDN);
74                 printf ("\n     got ");
75                 mpfr_out_str (stdout, 10, 0, r2, MPFR_RNDN);
76                 printf ("\n");
77                 exit (1);
78               }
79           }
80 
81   mpfr_clears (a, b, c, r1, r2, (mpfr_ptr) 0);
82 }
83 
84 static void
85 test_overflow1 (void)
86 {
87   mpfr_t x, y, z, r;
88   int inex;
89 
90   mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0);
91   MPFR_SET_POS (x);
92   mpfr_setmax (x, mpfr_get_emax ());  /* x = 2^emax - ulp */
93   mpfr_set_ui (y, 2, MPFR_RNDN);       /* y = 2 */
94   mpfr_neg (z, x, MPFR_RNDN);          /* z = -x = -(2^emax - ulp) */
95   mpfr_clear_flags ();
96   /* The intermediate multiplication x * y overflows, but x * y + z = x
97      is representable. */
98   inex = mpfr_fma (r, x, y, z, MPFR_RNDN);
99   if (inex || ! mpfr_equal_p (r, x))
100     {
101       printf ("Error in test_overflow1\nexpected ");
102       mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
103       printf (" with inex = 0\n     got ");
104       mpfr_out_str (stdout, 2, 0, r, MPFR_RNDN);
105       printf (" with inex = %d\n", inex);
106       exit (1);
107     }
108   if (mpfr_overflow_p ())
109     {
110       printf ("Error in test_overflow1: overflow flag set\n");
111       exit (1);
112     }
113   mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
114 }
115 
116 static void
117 test_overflow2 (void)
118 {
119   mpfr_t x, y, z, r;
120   int i, inex, rnd, err = 0;
121 
122   mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0);
123 
124   MPFR_SET_POS (x);
125   mpfr_setmin (x, mpfr_get_emax ());  /* x = 0.1@emax */
126   mpfr_set_si (y, -2, MPFR_RNDN);      /* y = -2 */
127   /* The intermediate multiplication x * y will overflow. */
128 
129   for (i = -9; i <= 9; i++)
130     RND_LOOP (rnd)
131       {
132         int inf, overflow;
133 
134         inf = rnd == MPFR_RNDN || rnd == MPFR_RNDD || rnd == MPFR_RNDA;
135         overflow = inf || i <= 0;
136 
137         inex = mpfr_set_si_2exp (z, i, mpfr_get_emin (), MPFR_RNDN);
138         MPFR_ASSERTN (inex == 0);
139 
140         mpfr_clear_flags ();
141         /* One has: x * y = -1@emax exactly (but not representable). */
142         inex = mpfr_fma (r, x, y, z, (mpfr_rnd_t) rnd);
143         if (overflow ^ (mpfr_overflow_p () != 0))
144           {
145             printf ("Error in test_overflow2 (i = %d, %s): wrong overflow"
146                     " flag (should be %d)\n", i,
147                     mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), overflow);
148             err = 1;
149           }
150         if (mpfr_nanflag_p ())
151           {
152             printf ("Error in test_overflow2 (i = %d, %s): NaN flag should"
153                     " not be set\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
154             err = 1;
155           }
156         if (mpfr_nan_p (r))
157           {
158             printf ("Error in test_overflow2 (i = %d, %s): got NaN\n",
159                     i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
160             err = 1;
161           }
162         else if (MPFR_SIGN (r) >= 0)
163           {
164             printf ("Error in test_overflow2 (i = %d, %s): wrong sign "
165                     "(+ instead of -)\n", i,
166                     mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
167             err = 1;
168           }
169         else if (inf && ! mpfr_inf_p (r))
170           {
171             printf ("Error in test_overflow2 (i = %d, %s): expected -Inf,"
172                     " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
173             mpfr_dump (r);
174             err = 1;
175           }
176         else if (!inf && (mpfr_inf_p (r) ||
177                           (mpfr_nextbelow (r), ! mpfr_inf_p (r))))
178           {
179             printf ("Error in test_overflow2 (i = %d, %s): expected -MAX,"
180                     " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
181             mpfr_dump (r);
182             err = 1;
183           }
184         if (inf ? inex >= 0 : inex <= 0)
185           {
186             printf ("Error in test_overflow2 (i = %d, %s): wrong inexact"
187                     " flag (got %d)\n", i,
188                     mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex);
189             err = 1;
190           }
191 
192       }
193 
194   if (err)
195     exit (1);
196   mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
197 }
198 
199 static void
200 test_underflow1 (void)
201 {
202   mpfr_t x, y, z, r;
203   int inex, signy, signz, rnd, err = 0;
204 
205   mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0);
206 
207   MPFR_SET_POS (x);
208   mpfr_setmin (x, mpfr_get_emin ());  /* x = 0.1@emin */
209 
210   for (signy = -1; signy <= 1; signy += 2)
211     {
212       mpfr_set_si_2exp (y, signy, -1, MPFR_RNDN);  /* |y| = 1/2 */
213       for (signz = -3; signz <= 3; signz += 2)
214         {
215           RND_LOOP (rnd)
216             {
217               mpfr_set_si (z, signz, MPFR_RNDN);
218               if (ABS (signz) != 1)
219                 mpfr_setmax (z, mpfr_get_emax ());
220               /* |z| = 1 or 2^emax - ulp */
221               mpfr_clear_flags ();
222               inex = mpfr_fma (r, x, y, z, (mpfr_rnd_t) rnd);
223 #define ERRTU1 "Error in test_underflow1 (signy = %d, signz = %d, %s)\n  "
224               if (mpfr_nanflag_p ())
225                 {
226                   printf (ERRTU1 "NaN flag is set\n", signy, signz,
227                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
228                   err = 1;
229                 }
230               if (signy < 0 && MPFR_IS_LIKE_RNDD(rnd, signz))
231                 mpfr_nextbelow (z);
232               if (signy > 0 && MPFR_IS_LIKE_RNDU(rnd, signz))
233                 mpfr_nextabove (z);
234               if ((mpfr_overflow_p () != 0) ^ (mpfr_inf_p (z) != 0))
235                 {
236                   printf (ERRTU1 "wrong overflow flag\n", signy, signz,
237                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
238                   err = 1;
239                 }
240               if (mpfr_underflow_p ())
241                 {
242                   printf (ERRTU1 "underflow flag is set\n", signy, signz,
243                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
244                   err = 1;
245                 }
246               if (! mpfr_equal_p (r, z))
247                 {
248                   printf (ERRTU1 "got ", signy, signz,
249                           mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
250                   mpfr_print_binary (r);
251                   printf (" instead of ");
252                   mpfr_print_binary (z);
253                   printf ("\n");
254                   err = 1;
255                 }
256               if (inex >= 0 && (rnd == MPFR_RNDD ||
257                                 (rnd == MPFR_RNDZ && signz > 0) ||
258                                 (rnd == MPFR_RNDN && signy > 0)))
259                 {
260                   printf (ERRTU1 "ternary value = %d instead of < 0\n",
261                           signy, signz, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
262                           inex);
263                   err = 1;
264                 }
265               if (inex <= 0 && (rnd == MPFR_RNDU ||
266                                 (rnd == MPFR_RNDZ && signz < 0) ||
267                                 (rnd == MPFR_RNDN && signy < 0)))
268                 {
269                   printf (ERRTU1 "ternary value = %d instead of > 0\n",
270                           signy, signz, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
271                           inex);
272                   err = 1;
273                 }
274             }
275         }
276     }
277 
278   if (err)
279     exit (1);
280   mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
281 }
282 
283 static void
284 test_underflow2 (void)
285 {
286   mpfr_t x, y, z, r;
287   int b, i, inex, same, err = 0;
288 
289   mpfr_inits2 (32, x, y, z, r, (mpfr_ptr) 0);
290 
291   mpfr_set_si_2exp (z, 1, mpfr_get_emin (), MPFR_RNDN);   /* z = 2^emin */
292   mpfr_set_si_2exp (x, 1, mpfr_get_emin (), MPFR_RNDN);   /* x = 2^emin */
293 
294   for (b = 0; b <= 1; b++)
295     {
296       for (i = 15; i <= 17; i++)
297         {
298           mpfr_set_si_2exp (y, i, -4 - MPFR_PREC (z), MPFR_RNDN);
299           /*  z = 1.000...00b
300            * xy =            01111
301            *   or            10000
302            *   or            10001
303            */
304           mpfr_clear_flags ();
305           inex = mpfr_fma (r, x, y, z, MPFR_RNDN);
306 #define ERRTU2 "Error in test_underflow2 (b = %d, i = %d)\n  "
307           if (__gmpfr_flags != MPFR_FLAGS_INEXACT)
308             {
309               printf (ERRTU2 "flags = %u instead of %u\n", b, i,
310                       __gmpfr_flags, (unsigned int) MPFR_FLAGS_INEXACT);
311               err = 1;
312             }
313           same = i == 15 || (i == 16 && b == 0);
314           if (same ? (inex >= 0) : (inex <= 0))
315             {
316               printf (ERRTU2 "incorrect ternary value (%d instead of %c 0)\n",
317                       b, i, inex, same ? '<' : '>');
318               err = 1;
319             }
320           mpfr_set (y, z, MPFR_RNDN);
321           if (!same)
322             mpfr_nextabove (y);
323           if (! mpfr_equal_p (r, y))
324             {
325               printf (ERRTU2 "expected ", b, i);
326               mpfr_dump (y);
327               printf ("  got      ");
328               mpfr_dump (r);
329               err = 1;
330             }
331         }
332       mpfr_nextabove (z);
333     }
334 
335   if (err)
336     exit (1);
337   mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
338 }
339 
340 static void
341 bug20101018 (void)
342 {
343   mpfr_t x, y, z, t, u;
344   int i;
345 
346   mpfr_init2 (x, 64);
347   mpfr_init2 (y, 64);
348   mpfr_init2 (z, 64);
349   mpfr_init2 (t, 64);
350   mpfr_init2 (u, 64);
351 
352   mpfr_set_str (x, "0xf.fffffffffffffffp-14766", 16, MPFR_RNDN);
353   mpfr_set_str (y, "-0xf.fffffffffffffffp+317", 16, MPFR_RNDN);
354   mpfr_set_str (z, "0x8.3ffffffffffe3ffp-14443", 16, MPFR_RNDN);
355   mpfr_set_str (t, "0x8.7ffffffffffc7ffp-14444", 16, MPFR_RNDN);
356   i = mpfr_fma (u, x, y, z, MPFR_RNDN);
357   if (mpfr_cmp (u, t) != 0)
358     {
359       printf ("Wrong result in bug20101018 (a)\n");
360       printf ("Expected ");
361       mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
362       printf ("\nGot      ");
363       mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
364       printf ("\n");
365       exit (1);
366     }
367   if (i <= 0)
368     {
369       printf ("Wrong ternary value in bug20101018 (a)\n");
370       printf ("Expected > 0\n");
371       printf ("Got      %d\n", i);
372       exit (1);
373     }
374 
375   mpfr_set_str (x, "-0xf.fffffffffffffffp-11420", 16, MPFR_RNDN);
376   mpfr_set_str (y, "0xf.fffffffffffffffp+9863", 16, MPFR_RNDN);
377   mpfr_set_str (z, "0x8.fffff80ffffffffp-1551", 16, MPFR_RNDN);
378   mpfr_set_str (t, "0x9.fffff01ffffffffp-1552", 16, MPFR_RNDN);
379   i = mpfr_fma (u, x, y, z, MPFR_RNDN);
380   if (mpfr_cmp (u, t) != 0)
381     {
382       printf ("Wrong result in bug20101018 (b)\n");
383       printf ("Expected ");
384       mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
385       printf ("\nGot      ");
386       mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
387       printf ("\n");
388       exit (1);
389     }
390   if (i <= 0)
391     {
392       printf ("Wrong ternary value in bug20101018 (b)\n");
393       printf ("Expected > 0\n");
394       printf ("Got      %d\n", i);
395       exit (1);
396     }
397 
398   mpfr_set_str (x, "0xf.fffffffffffffffp-2125", 16, MPFR_RNDN);
399   mpfr_set_str (y, "-0xf.fffffffffffffffp-6000", 16, MPFR_RNDN);
400   mpfr_set_str (z, "0x8p-8119", 16, MPFR_RNDN);
401   mpfr_set_str (t, "0x8.000000000000001p-8120", 16, MPFR_RNDN);
402   i = mpfr_fma (u, x, y, z, MPFR_RNDN);
403   if (mpfr_cmp (u, t) != 0)
404     {
405       printf ("Wrong result in bug20101018 (c)\n");
406       printf ("Expected ");
407       mpfr_out_str (stdout, 16, 0, t, MPFR_RNDN);
408       printf ("\nGot      ");
409       mpfr_out_str (stdout, 16, 0, u, MPFR_RNDN);
410       printf ("\n");
411       exit (1);
412     }
413   if (i <= 0)
414     {
415       printf ("Wrong ternary value in bug20101018 (c)\n");
416       printf ("Expected > 0\n");
417       printf ("Got      %d\n", i);
418       exit (1);
419     }
420 
421   mpfr_clear (x);
422   mpfr_clear (y);
423   mpfr_clear (z);
424   mpfr_clear (t);
425   mpfr_clear (u);
426 }
427 
428 int
429 main (int argc, char *argv[])
430 {
431   mpfr_t x, y, z, s;
432   MPFR_SAVE_EXPO_DECL (expo);
433 
434   tests_start_mpfr ();
435 
436   bug20101018 ();
437 
438   mpfr_init (x);
439   mpfr_init (s);
440   mpfr_init (y);
441   mpfr_init (z);
442 
443   /* check special cases */
444   mpfr_set_prec (x, 2);
445   mpfr_set_prec (y, 2);
446   mpfr_set_prec (z, 2);
447   mpfr_set_prec (s, 2);
448   mpfr_set_str (x, "-0.75", 10, MPFR_RNDN);
449   mpfr_set_str (y, "0.5", 10, MPFR_RNDN);
450   mpfr_set_str (z, "0.375", 10, MPFR_RNDN);
451   mpfr_fma (s, x, y, z, MPFR_RNDU); /* result is 0 */
452   if (mpfr_cmp_ui(s, 0))
453     {
454       printf("Error: -0.75 * 0.5 + 0.375 should be equal to 0 for prec=2\n");
455       exit(1);
456     }
457 
458   mpfr_set_prec (x, 27);
459   mpfr_set_prec (y, 27);
460   mpfr_set_prec (z, 27);
461   mpfr_set_prec (s, 27);
462   mpfr_set_str_binary (x, "1.11111111111111111111111111e-1");
463   mpfr_set (y, x, MPFR_RNDN);
464   mpfr_set_str_binary (z, "-1.00011110100011001011001001e-1");
465   if (mpfr_fma (s, x, y, z, MPFR_RNDN) >= 0)
466     {
467       printf ("Wrong inexact flag for x=y=1-2^(-27)\n");
468       exit (1);
469     }
470 
471   mpfr_set_nan (x);
472   mpfr_urandomb (y, RANDS);
473   mpfr_urandomb (z, RANDS);
474   mpfr_fma (s, x, y, z, MPFR_RNDN);
475   if (!mpfr_nan_p (s))
476     {
477       printf ("evaluation of function in x=NAN does not return NAN");
478       exit (1);
479     }
480 
481   mpfr_set_nan (y);
482   mpfr_urandomb (x, RANDS);
483   mpfr_urandomb (z, RANDS);
484   mpfr_fma (s, x, y, z, MPFR_RNDN);
485   if (!mpfr_nan_p(s))
486     {
487       printf ("evaluation of function in y=NAN does not return NAN");
488       exit (1);
489     }
490 
491   mpfr_set_nan (z);
492   mpfr_urandomb (y, RANDS);
493   mpfr_urandomb (x, RANDS);
494   mpfr_fma (s, x, y, z, MPFR_RNDN);
495   if (!mpfr_nan_p (s))
496     {
497       printf ("evaluation of function in z=NAN does not return NAN");
498       exit (1);
499     }
500 
501   mpfr_set_inf (x, 1);
502   mpfr_set_inf (y, 1);
503   mpfr_set_inf (z, 1);
504   mpfr_fma (s, x, y, z, MPFR_RNDN);
505   if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
506     {
507       printf ("Error for (+inf) * (+inf) + (+inf)\n");
508       exit (1);
509     }
510 
511   mpfr_set_inf (x, -1);
512   mpfr_set_inf (y, -1);
513   mpfr_set_inf (z, 1);
514   mpfr_fma (s, x, y, z, MPFR_RNDN);
515   if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
516     {
517       printf ("Error for (-inf) * (-inf) + (+inf)\n");
518       exit (1);
519     }
520 
521   mpfr_set_inf (x, 1);
522   mpfr_set_inf (y, -1);
523   mpfr_set_inf (z, -1);
524   mpfr_fma (s, x, y, z, MPFR_RNDN);
525   if (!mpfr_inf_p (s) || mpfr_sgn (s) > 0)
526     {
527       printf ("Error for (+inf) * (-inf) + (-inf)\n");
528       exit (1);
529     }
530 
531   mpfr_set_inf (x, -1);
532   mpfr_set_inf (y, 1);
533   mpfr_set_inf (z, -1);
534   mpfr_fma (s, x, y, z, MPFR_RNDN);
535   if (!mpfr_inf_p (s) || mpfr_sgn (s) > 0)
536     {
537       printf ("Error for (-inf) * (+inf) + (-inf)\n");
538       exit (1);
539     }
540 
541   mpfr_set_inf (x, 1);
542   mpfr_set_ui (y, 0, MPFR_RNDN);
543   mpfr_urandomb (z, RANDS);
544   mpfr_fma (s, x, y, z, MPFR_RNDN);
545   if (!mpfr_nan_p (s))
546     {
547       printf ("evaluation of function in x=INF y=0  does not return NAN");
548       exit (1);
549     }
550 
551   mpfr_set_inf (y, 1);
552   mpfr_set_ui (x, 0, MPFR_RNDN);
553   mpfr_urandomb (z, RANDS);
554   mpfr_fma (s, x, y, z, MPFR_RNDN);
555   if (!mpfr_nan_p (s))
556     {
557       printf ("evaluation of function in x=0 y=INF does not return NAN");
558       exit (1);
559     }
560 
561   mpfr_set_inf (x, 1);
562   mpfr_urandomb (y, RANDS); /* always positive */
563   mpfr_set_inf (z, -1);
564   mpfr_fma (s, x, y, z, MPFR_RNDN);
565   if (!mpfr_nan_p (s))
566     {
567       printf ("evaluation of function in x=INF y>0 z=-INF does not return NAN");
568       exit (1);
569     }
570 
571   mpfr_set_inf (y, 1);
572   mpfr_urandomb (x, RANDS);
573   mpfr_set_inf (z, -1);
574   mpfr_fma (s, x, y, z, MPFR_RNDN);
575   if (!mpfr_nan_p (s))
576     {
577       printf ("evaluation of function in x>0 y=INF z=-INF does not return NAN");
578       exit (1);
579     }
580 
581   mpfr_set_inf (x, 1);
582   mpfr_urandomb (y, RANDS);
583   mpfr_urandomb (z, RANDS);
584   mpfr_fma (s, x, y, z, MPFR_RNDN);
585   if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
586     {
587       printf ("evaluation of function in x=INF does not return INF");
588       exit (1);
589     }
590 
591   mpfr_set_inf (y, 1);
592   mpfr_urandomb (x, RANDS);
593   mpfr_urandomb (z, RANDS);
594   mpfr_fma (s, x, y, z, MPFR_RNDN);
595   if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
596     {
597       printf ("evaluation of function in y=INF does not return INF");
598       exit (1);
599     }
600 
601   mpfr_set_inf (z, 1);
602   mpfr_urandomb (x, RANDS);
603   mpfr_urandomb (y, RANDS);
604   mpfr_fma (s, x, y, z, MPFR_RNDN);
605   if (!mpfr_inf_p (s) || mpfr_sgn (s) < 0)
606     {
607       printf ("evaluation of function in z=INF does not return INF");
608       exit (1);
609     }
610 
611   mpfr_set_ui (x, 0, MPFR_RNDN);
612   mpfr_urandomb (y, RANDS);
613   mpfr_urandomb (z, RANDS);
614   mpfr_fma (s, x, y, z, MPFR_RNDN);
615   if (mpfr_cmp (s, z))
616     {
617       printf ("evaluation of function in x=0 does not return z\n");
618       exit (1);
619     }
620 
621   mpfr_set_ui (y, 0, MPFR_RNDN);
622   mpfr_urandomb (x, RANDS);
623   mpfr_urandomb (z, RANDS);
624   mpfr_fma (s, x, y, z, MPFR_RNDN);
625   if (mpfr_cmp (s, z))
626     {
627       printf ("evaluation of function in y=0 does not return z\n");
628       exit (1);
629     }
630 
631   {
632     mpfr_prec_t prec;
633     mpfr_t t, slong;
634     mpfr_rnd_t rnd;
635     int inexact, compare;
636     unsigned int n;
637 
638     mpfr_prec_t p0=2, p1=200;
639     unsigned int N=200;
640 
641     mpfr_init (t);
642     mpfr_init (slong);
643 
644     /* generic test */
645     for (prec = p0; prec <= p1; prec++)
646     {
647       mpfr_set_prec (x, prec);
648       mpfr_set_prec (y, prec);
649       mpfr_set_prec (z, prec);
650       mpfr_set_prec (s, prec);
651       mpfr_set_prec (t, prec);
652 
653       for (n=0; n<N; n++)
654         {
655           mpfr_urandomb (x, RANDS);
656           mpfr_urandomb (y, RANDS);
657           mpfr_urandomb (z, RANDS);
658 
659           if (randlimb () % 2)
660             mpfr_neg (x, x, MPFR_RNDN);
661           if (randlimb () % 2)
662             mpfr_neg (y, y, MPFR_RNDN);
663           if (randlimb () % 2)
664             mpfr_neg (z, z, MPFR_RNDN);
665 
666           rnd = RND_RAND ();
667           mpfr_set_prec (slong, 2 * prec);
668           if (mpfr_mul (slong, x, y, rnd))
669             {
670               printf ("x*y should be exact\n");
671               exit (1);
672             }
673           compare = mpfr_add (t, slong, z, rnd);
674           inexact = mpfr_fma (s, x, y, z, rnd);
675           if (mpfr_cmp (s, t))
676             {
677               printf ("results differ for x=");
678               mpfr_out_str (stdout, 2, prec, x, MPFR_RNDN);
679               printf ("  y=");
680               mpfr_out_str (stdout, 2, prec, y, MPFR_RNDN);
681               printf ("  z=");
682               mpfr_out_str (stdout, 2, prec, z, MPFR_RNDN);
683               printf (" prec=%u rnd_mode=%s\n", (unsigned int) prec,
684                       mpfr_print_rnd_mode (rnd));
685               printf ("got      ");
686               mpfr_out_str (stdout, 2, prec, s, MPFR_RNDN);
687               puts ("");
688               printf ("expected ");
689               mpfr_out_str (stdout, 2, prec, t, MPFR_RNDN);
690               puts ("");
691               printf ("approx  ");
692               mpfr_print_binary (slong);
693               puts ("");
694               exit (1);
695             }
696           if (((inexact == 0) && (compare != 0)) ||
697               ((inexact < 0) && (compare >= 0)) ||
698               ((inexact > 0) && (compare <= 0)))
699             {
700               printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n",
701                       mpfr_print_rnd_mode (rnd), compare, inexact);
702               printf (" x="); mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
703               printf (" y="); mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
704               printf (" z="); mpfr_out_str (stdout, 2, 0, z, MPFR_RNDN);
705               printf (" s="); mpfr_out_str (stdout, 2, 0, s, MPFR_RNDN);
706               printf ("\n");
707               exit (1);
708             }
709         }
710     }
711   mpfr_clear (t);
712   mpfr_clear (slong);
713 
714   }
715   mpfr_clear (x);
716   mpfr_clear (y);
717   mpfr_clear (z);
718   mpfr_clear (s);
719 
720   test_exact ();
721 
722   MPFR_SAVE_EXPO_MARK (expo);
723   test_overflow1 ();
724   test_overflow2 ();
725   test_underflow1 ();
726   test_underflow2 ();
727   MPFR_SAVE_EXPO_FREE (expo);
728 
729   tests_end_mpfr ();
730   return 0;
731 }
732