← Back to branch summary

~ubuntu-branches/debian/sid/genius/sid

« back to all changes in this revision

Viewing changes to mpfr/tests/trint.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel Holbach
  • Date: 2006-08-21 12:57:45 UTC
  • Revision ID: james.westby@ubuntu.com-20060821125745-sl9ks8v7fq324bdf
Tags: upstream-0.7.6.1
Import upstream version 0.7.6.1

Show diffs side-by-side

added added

removed removed

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