← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/mpfr/tests/tsub_ui.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2007-04-09 11:51:51 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20070409115151-ql8cr0kalzx1jmla
Tags: 0.9i-20070324-2
Upload to unstable. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpfr/tests/tsub_ui.c
1
 
/* Test file for mpfr_sub_ui
2
 
 
3
 
Copyright 2000, 2001, 2002 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., 59 Temple Place - Suite 330, Boston,
20
 
MA 02111-1307, USA. */
21
 
 
22
 
#include <math.h>
23
 
#include <stdio.h>
24
 
#include <stdlib.h>
25
 
#include <time.h>
26
 
#include "gmp.h"
27
 
#include "mpfr.h"
28
 
#include "mpfr-impl.h"
29
 
#include "mpfr-test.h"
30
 
 
31
 
void check_two_sum _PROTO ((mp_prec_t));
32
 
void check3 _PROTO ((double, unsigned long, mp_rnd_t, double));
33
 
 
34
 
#define check(x,y,r) check3(x,y,r,0.0)
35
 
 
36
 
/* checks that x+y gives the same results in double
37
 
   and with mpfr with 53 bits of precision */
38
 
void
39
 
check3 (double x, unsigned long y, mp_rnd_t rnd_mode, double z1)
40
 
{
41
 
  double z2;
42
 
  mpfr_t xx,zz;
43
 
 
44
 
  mpfr_init(xx);
45
 
  mpfr_init(zz);
46
 
  mpfr_set_prec(xx, 53);
47
 
  mpfr_set_prec(zz, 53);
48
 
  mpfr_set_d(xx, x, rnd_mode);
49
 
  mpfr_sub_ui(zz, xx, y, rnd_mode);
50
 
#ifdef MPFR_HAVE_FESETROUND
51
 
  mpfr_set_machine_rnd_mode(rnd_mode);
52
 
#endif
53
 
  if (z1==0.0) z1 = x-y;
54
 
  z2 = mpfr_get_d1 (zz);
55
 
  if (z1!=z2 && !(isnan(z1) && isnan(z2))) {
56
 
    printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
57
 
    printf("mpfr_sub_ui failed for x=%1.20e y=%lu with rnd_mode=%s\n",
58
 
           x, y, mpfr_print_rnd_mode(rnd_mode));
59
 
    exit(1);
60
 
  }
61
 
  mpfr_clear(xx);
62
 
  mpfr_clear(zz);
63
 
}
64
 
 
65
 
/* FastTwoSum: if EXP(x) >= EXP(y), u = o(x+y), v = o(u-x), w = o(y-v),
66
 
               then x + y = u + w
67
 
thus if u = o(y-x), v = o(u+x), w = o(v-y), then y-x = u-w */
68
 
void
69
 
check_two_sum (mp_prec_t p)
70
 
{
71
 
  unsigned int x;
72
 
  mpfr_t y, u, v, w;
73
 
  mp_rnd_t rnd;
74
 
  int inexact;
75
 
  
76
 
  mpfr_init2 (y, p);
77
 
  mpfr_init2 (u, p);
78
 
  mpfr_init2 (v, p);
79
 
  mpfr_init2 (w, p);
80
 
  do
81
 
    {
82
 
      x = LONG_RAND ();
83
 
    }
84
 
  while (x < 1);
85
 
  mpfr_random (y);
86
 
  rnd = LONG_RAND() % 4;
87
 
  rnd = GMP_RNDN;
88
 
  inexact = mpfr_sub_ui (u, y, x, GMP_RNDN);
89
 
  mpfr_add_ui (v, u, x, GMP_RNDN);
90
 
  mpfr_sub (w, v, y, GMP_RNDN);
91
 
  /* as u - (y-x) = w, we should have inexact and w of same sign */
92
 
  if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
93
 
      ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) ||
94
 
      ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0)))
95
 
    {
96
 
      fprintf (stderr, "Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p,
97
 
               mpfr_print_rnd_mode (rnd));
98
 
      printf ("x=%u\n", x);
99
 
      printf ("y="); mpfr_print_binary(y); putchar('\n');
100
 
      printf ("u="); mpfr_print_binary(u); putchar('\n');
101
 
      printf ("v="); mpfr_print_binary(v); putchar('\n');
102
 
      printf ("w="); mpfr_print_binary(w); putchar('\n');
103
 
      printf ("inexact = %d\n", inexact);
104
 
      exit (1);
105
 
    }
106
 
  mpfr_clear (y);
107
 
  mpfr_clear (u);
108
 
  mpfr_clear (v);
109
 
  mpfr_clear (w);
110
 
}
111
 
 
112
 
int
113
 
main (int argc, char *argv[])
114
 
{
115
 
  mp_prec_t p;
116
 
  int k;
117
 
#ifdef MPFR_HAVE_FESETROUND
118
 
  double x; unsigned long y, N; int i,rnd_mode,rnd;
119
 
 
120
 
  mpfr_test_init ();
121
 
 
122
 
  SEED_RAND (time(NULL));
123
 
  N = (argc<2) ? 1000000 : atoi(argv[1]);
124
 
  rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
125
 
  for (i=0;i<1000000;i++) {
126
 
    x = drand();
127
 
    y = LONG_RAND();
128
 
    if (ABS(x)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
129
 
      /* avoid denormalized numbers and overflows */
130
 
      rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
131
 
      check(x, y, rnd);
132
 
    }
133
 
  } 
134
 
#endif
135
 
  
136
 
  for (p=2; p<200; p++)
137
 
    for (k=0; k<200; k++)
138
 
      check_two_sum (p);
139
 
 
140
 
  check3 (0.9999999999, 1, GMP_RNDN, -56295.0 / 562949953421312.0);
141
 
#ifdef HAVE_INFS
142
 
  check3 (DBL_NAN, 1, GMP_RNDN, DBL_NAN);
143
 
  check3 (DBL_POS_INF, 1, GMP_RNDN, DBL_POS_INF);
144
 
  check3 (DBL_NEG_INF, 1, GMP_RNDN, DBL_NEG_INF);
145
 
#endif
146
 
 
147
 
  return 0;
148
 
}
149
 
 
150