1/* Test file for mpfr_rint, mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round.
2
3Copyright 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 <stdlib.h>
24
25#include "mpfr-test.h"
26
27#if __MPFR_STDC (199901L)
28# include <math.h>
29#endif
30
31static void
32special (void)
33{
34  mpfr_t x, y;
35  mpfr_exp_t emax;
36
37  mpfr_init (x);
38  mpfr_init (y);
39
40  mpfr_set_nan (x);
41  mpfr_rint (y, x, MPFR_RNDN);
42  MPFR_ASSERTN(mpfr_nan_p (y));
43
44  mpfr_set_inf (x, 1);
45  mpfr_rint (y, x, MPFR_RNDN);
46  MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0);
47
48  mpfr_set_inf (x, -1);
49  mpfr_rint (y, x, MPFR_RNDN);
50  MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) < 0);
51
52  mpfr_set_ui (x, 0, MPFR_RNDN);
53  mpfr_rint (y, x, MPFR_RNDN);
54  MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y));
55
56  mpfr_set_ui (x, 0, MPFR_RNDN);
57  mpfr_neg (x, x, MPFR_RNDN);
58  mpfr_rint (y, x, MPFR_RNDN);
59  MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_NEG(y));
60
61  /* coverage test */
62  mpfr_set_prec (x, 2);
63  mpfr_set_ui (x, 1, MPFR_RNDN);
64  mpfr_mul_2exp (x, x, mp_bits_per_limb, MPFR_RNDN);
65  mpfr_rint (y, x, MPFR_RNDN);
66  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
67
68  /* another coverage test */
69  emax = mpfr_get_emax ();
70  set_emax (1);
71  mpfr_set_prec (x, 3);
72  mpfr_set_str_binary (x, "1.11E0");
73  mpfr_set_prec (y, 2);
74  mpfr_rint (y, x, MPFR_RNDU); /* x rounds to 1.0E1=0.1E2 which overflows */
75  MPFR_ASSERTN(mpfr_inf_p (y) && mpfr_sgn (y) > 0);
76  set_emax (emax);
77
78  /* yet another */
79  mpfr_set_prec (x, 97);
80  mpfr_set_prec (y, 96);
81  mpfr_set_str_binary (x, "-0.1011111001101111000111011100011100000110110110110000000111010001000101001111101010101011010111100E97");
82  mpfr_rint (y, x, MPFR_RNDN);
83  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
84
85  mpfr_set_prec (x, 53);
86  mpfr_set_prec (y, 53);
87  mpfr_set_str_binary (x, "0.10101100000000101001010101111111000000011111010000010E-1");
88  mpfr_rint (y, x, MPFR_RNDU);
89  MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0);
90  mpfr_rint (y, x, MPFR_RNDD);
91  MPFR_ASSERTN(mpfr_cmp_ui (y, 0) == 0 && MPFR_IS_POS(y));
92
93  mpfr_set_prec (x, 36);
94  mpfr_set_prec (y, 2);
95  mpfr_set_str_binary (x, "-11000110101010111111110111001.0000100");
96  mpfr_rint (y, x, MPFR_RNDN);
97  mpfr_set_str_binary (x, "-11E27");
98  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
99
100  mpfr_set_prec (x, 39);
101  mpfr_set_prec (y, 29);
102  mpfr_set_str_binary (x, "-0.100010110100011010001111001001001100111E39");
103  mpfr_rint (y, x, MPFR_RNDN);
104  mpfr_set_str_binary (x, "-0.10001011010001101000111100101E39");
105  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
106
107  mpfr_set_prec (x, 46);
108  mpfr_set_prec (y, 32);
109  mpfr_set_str_binary (x, "-0.1011100110100101000001011111101011001001101001E32");
110  mpfr_rint (y, x, MPFR_RNDN);
111  mpfr_set_str_binary (x, "-0.10111001101001010000010111111011E32");
112  MPFR_ASSERTN(mpfr_cmp (y, x) == 0);
113
114  /* coverage test for mpfr_round */
115  mpfr_set_prec (x, 3);
116  mpfr_set_str_binary (x, "1.01E1"); /* 2.5 */
117  mpfr_set_prec (y, 2);
118  mpfr_round (y, x);
119  /* since mpfr_round breaks ties away, should give 3 and not 2 as with
120     the "round to even" rule */
121  MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0);
122  /* same test for the function */
123  (mpfr_round) (y, x);
124  MPFR_ASSERTN(mpfr_cmp_ui (y, 3) == 0);
125
126  mpfr_set_prec (x, 6);
127  mpfr_set_prec (y, 3);
128  mpfr_set_str_binary (x, "110.111");
129  mpfr_round (y, x);
130  if (mpfr_cmp_ui (y, 7))
131    {
132      printf ("Error in round(110.111)\n");
133      exit (1);
134    }
135
136  /* Bug found by  Mark J Watkins */
137  mpfr_set_prec (x, 84);
138  mpfr_set_str_binary (x,
139   "0.110011010010001000000111101101001111111100101110010000000000000" \
140                       "000000000000000000000E32");
141  mpfr_round (x, x);
142  if (mpfr_cmp_str (x, "0.1100110100100010000001111011010100000000000000" \
143                    "00000000000000000000000000000000000000E32", 2, MPFR_RNDN))
144    {
145      printf ("Rounding error when dest=src\n");
146      exit (1);
147    }
148
149  mpfr_clear (x);
150  mpfr_clear (y);
151}
152
153#if __MPFR_STDC (199901L)
154
155static void
156test_fct (double (*f)(double), int (*g)(), char *s, mpfr_rnd_t r)
157{
158  double d, y;
159  mpfr_t dd, yy;
160
161  mpfr_init2 (dd, 53);
162  mpfr_init2 (yy, 53);
163  for (d = -5.0; d <= 5.0; d += 0.25)
164    {
165      mpfr_set_d (dd, d, r);
166      y = (*f)(d);
167      if (g == &mpfr_rint)
168        mpfr_rint (yy, dd, r);
169      else
170        (*g)(yy, dd);
171      mpfr_set_d (dd, y, r);
172      if (mpfr_cmp (yy, dd))
173        {
174          printf ("test_against_libc: incorrect result for %s, rnd = %s,"
175                  " d = %g\ngot ", s, mpfr_print_rnd_mode (r), d);
176          mpfr_out_str (stdout, 10, 0, yy, MPFR_RNDN);
177          printf (" instead of %g\n", y);
178          exit (1);
179        }
180    }
181  mpfr_clear (dd);
182  mpfr_clear (yy);
183}
184
185#define TEST_FCT(F) test_fct (&F, &mpfr_##F, #F, r)
186
187static void
188test_against_libc (void)
189{
190  mpfr_rnd_t r = MPFR_RNDN;
191
192#if HAVE_ROUND
193  TEST_FCT (round);
194#endif
195#if HAVE_TRUNC
196  TEST_FCT (trunc);
197#endif
198#if HAVE_FLOOR
199  TEST_FCT (floor);
200#endif
201#if HAVE_CEIL
202  TEST_FCT (ceil);
203#endif
204#if HAVE_NEARBYINT
205  for (r = 0; r < MPFR_RND_MAX ; r++)
206    if (mpfr_set_machine_rnd_mode (r) == 0)
207      test_fct (&nearbyint, &mpfr_rint, "rint", r);
208#endif
209}
210
211#endif
212
213static void
214err (char *str, mp_size_t s, mpfr_t x, mpfr_t y, mpfr_prec_t p, mpfr_rnd_t r,
215     int trint, int inexact)
216{
217  printf ("Error: %s\ns = %u, p = %u, r = %s, trint = %d, inexact = %d\nx = ",
218          str, (unsigned int) s, (unsigned int) p, mpfr_print_rnd_mode (r),
219          trint, inexact);
220  mpfr_print_binary (x);
221  printf ("\ny = ");
222  mpfr_print_binary (y);
223  printf ("\n");
224  exit (1);
225}
226
227int
228main (int argc, char *argv[])
229{
230  mp_size_t s;
231  mpz_t z;
232  mpfr_prec_t p;
233  mpfr_t x, y, t, u, v;
234  int r;
235  int inexact, sign_t;
236
237  tests_start_mpfr ();
238
239  mpfr_init (x);
240  mpfr_init (y);
241  mpz_init (z);
242  mpfr_init (t);
243  mpfr_init (u);
244  mpfr_init (v);
245  mpz_set_ui (z, 1);
246  for (s = 2; s < 100; s++)
247    {
248      /* z has exactly s bits */
249
250      mpz_mul_2exp (z, z, 1);
251      if (randlimb () % 2)
252        mpz_add_ui (z, z, 1);
253      mpfr_set_prec (x, s);
254      mpfr_set_prec (t, s);
255      mpfr_set_prec (u, s);
256      if (mpfr_set_z (x, z, MPFR_RNDN))
257        {
258          printf ("Error: mpfr_set_z should be exact (s = %u)\n",
259                  (unsigned int) s);
260          exit (1);
261        }
262      if (randlimb () % 2)
263        mpfr_neg (x, x, MPFR_RNDN);
264      if (randlimb () % 2)
265        mpfr_div_2ui (x, x, randlimb () % s, MPFR_RNDN);
266      for (p = 2; p < 100; p++)
267        {
268          int trint;
269          mpfr_set_prec (y, p);
270          mpfr_set_prec (v, p);
271          for (r = 0; r < MPFR_RND_MAX ; r++)
272            for (trint = 0; trint < 3; trint++)
273              {
274                if (trint == 2)
275                  inexact = mpfr_rint (y, x, (mpfr_rnd_t) r);
276                else if (r == MPFR_RNDN)
277                  inexact = mpfr_round (y, x);
278                else if (r == MPFR_RNDZ)
279                  inexact = (trint ? mpfr_trunc (y, x) :
280                             mpfr_rint_trunc (y, x, MPFR_RNDZ));
281                else if (r == MPFR_RNDU)
282                  inexact = (trint ? mpfr_ceil (y, x) :
283                             mpfr_rint_ceil (y, x, MPFR_RNDU));
284                else /* r = MPFR_RNDD */
285                  inexact = (trint ? mpfr_floor (y, x) :
286                             mpfr_rint_floor (y, x, MPFR_RNDD));
287                if (mpfr_sub (t, y, x, MPFR_RNDN))
288                  err ("subtraction 1 should be exact",
289                       s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
290                sign_t = mpfr_cmp_ui (t, 0);
291                if (trint != 0 &&
292                    (((inexact == 0) && (sign_t != 0)) ||
293                     ((inexact < 0) && (sign_t >= 0)) ||
294                     ((inexact > 0) && (sign_t <= 0))))
295                  err ("wrong inexact flag", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
296                if (inexact == 0)
297                  continue; /* end of the test for exact results */
298
299                if (((r == MPFR_RNDD || (r == MPFR_RNDZ && MPFR_SIGN (x) > 0))
300                     && inexact > 0) ||
301                    ((r == MPFR_RNDU || (r == MPFR_RNDZ && MPFR_SIGN (x) < 0))
302                     && inexact < 0))
303                  err ("wrong rounding direction",
304                       s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
305                if (inexact < 0)
306                  {
307                    mpfr_add_ui (v, y, 1, MPFR_RNDU);
308                    if (mpfr_cmp (v, x) <= 0)
309                      err ("representable integer between x and its "
310                           "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
311                  }
312                else
313                  {
314                    mpfr_sub_ui (v, y, 1, MPFR_RNDD);
315                    if (mpfr_cmp (v, x) >= 0)
316                      err ("representable integer between x and its "
317                           "rounded value", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
318                  }
319                if (r == MPFR_RNDN)
320                  {
321                    int cmp;
322                    if (mpfr_sub (u, v, x, MPFR_RNDN))
323                      err ("subtraction 2 should be exact",
324                           s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
325                    cmp = mpfr_cmp_abs (t, u);
326                    if (cmp > 0)
327                      err ("faithful rounding, but not the nearest integer",
328                           s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
329                    if (cmp < 0)
330                      continue;
331                    /* |t| = |u|: x is the middle of two consecutive
332                       representable integers. */
333                    if (trint == 2)
334                      {
335                        /* halfway case for mpfr_rint in MPFR_RNDN rounding
336                           mode: round to an even integer or mantissa. */
337                        mpfr_div_2ui (y, y, 1, MPFR_RNDZ);
338                        if (!mpfr_integer_p (y))
339                          err ("halfway case for mpfr_rint, result isn't an"
340                               " even integer", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
341                        /* If floor(x) and ceil(x) aren't both representable
342                           integers, the mantissa must be even. */
343                        mpfr_sub (v, v, y, MPFR_RNDN);
344                        mpfr_abs (v, v, MPFR_RNDN);
345                        if (mpfr_cmp_ui (v, 1) != 0)
346                          {
347                            mpfr_div_2si (y, y, MPFR_EXP (y) - MPFR_PREC (y)
348                                          + 1, MPFR_RNDN);
349                            if (!mpfr_integer_p (y))
350                              err ("halfway case for mpfr_rint, mantissa isn't"
351                                   " even", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
352                          }
353                      }
354                    else
355                      { /* halfway case for mpfr_round: x must have been
356                           rounded away from zero. */
357                        if ((MPFR_SIGN (x) > 0 && inexact < 0) ||
358                            (MPFR_SIGN (x) < 0 && inexact > 0))
359                          err ("halfway case for mpfr_round, bad rounding"
360                               " direction", s, x, y, p, (mpfr_rnd_t) r, trint, inexact);
361                      }
362                  }
363              }
364        }
365    }
366  mpfr_clear (x);
367  mpfr_clear (y);
368  mpz_clear (z);
369  mpfr_clear (t);
370  mpfr_clear (u);
371  mpfr_clear (v);
372
373  special ();
374
375#if __MPFR_STDC (199901L)
376  if (argc > 1 && strcmp (argv[1], "-s") == 0)
377    test_against_libc ();
378#endif
379
380  tests_end_mpfr ();
381  return 0;
382}
383