1/* Test file for mpfr_fma.
2
3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by the Arenaire and Cacao projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 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. */
29static void
30test_exact (void)
31{
32  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
84static void
85test_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
116static void
117test_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
199static void
200test_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
283static void
284test_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
340static void
341bug20101018 (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
428int
429main (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