← 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.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.c
1
 
/* Test file for mpfr_sub.
2
 
 
3
 
Copyright 2001 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 "gmp.h"
26
 
#include "mpfr.h"
27
 
#include "mpfr-impl.h"
28
 
#include "mpfr-test.h"
29
 
 
30
 
void check_diverse _PROTO((void));
31
 
void bug_ddefour _PROTO((void));
32
 
void check_two_sum _PROTO((mp_prec_t));
33
 
void check_inexact _PROTO((void));
34
 
 
35
 
void
36
 
check_diverse (void)
37
 
{
38
 
  mpfr_t x, y, z;
39
 
  double res, got;
40
 
  int inexact;
41
 
 
42
 
  mpfr_init (x);
43
 
  mpfr_init (y);
44
 
  mpfr_init (z);
45
 
 
46
 
  mpfr_set_prec (x, 32);
47
 
  mpfr_set_prec (y, 63);
48
 
  mpfr_set_prec (z, 63);
49
 
  mpfr_set_str_raw (x, "0.101101111011011100100100100111E31");
50
 
  mpfr_set_str_raw (y, "0.111110010010100100110101101010001001100101110001000101110111111E-1");
51
 
  mpfr_sub (z, x, y, GMP_RNDN);
52
 
  mpfr_set_str_raw (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
53
 
  if (mpfr_cmp (z, y))
54
 
    {
55
 
      fprintf (stderr, "Error in mpfr_sub (5)\n");
56
 
      printf ("expected "); mpfr_print_binary (y); putchar ('\n');
57
 
      printf ("got      "); mpfr_print_binary (z); putchar ('\n');
58
 
      exit (1);
59
 
    }
60
 
 
61
 
  mpfr_set_prec (y, 63);
62
 
  mpfr_set_prec (z, 63);
63
 
  mpfr_set_str_raw (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
64
 
  mpfr_sub_ui (z, y, 1541116494, GMP_RNDN);
65
 
  mpfr_set_str_raw (y, "-0.11111001001010010011010110101E-1");
66
 
  if (mpfr_cmp (z, y))
67
 
    {
68
 
      fprintf (stderr, "Error in mpfr_sub (7)\n");
69
 
      printf ("expected "); mpfr_print_binary (y); putchar ('\n');
70
 
      printf ("got      "); mpfr_print_binary (z); putchar ('\n');
71
 
      exit (1);
72
 
    }
73
 
 
74
 
  mpfr_set_prec (y, 63);
75
 
  mpfr_set_prec (z, 63);
76
 
  mpfr_set_str_raw (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
77
 
  mpfr_sub_ui (z, y, 1541116494, GMP_RNDN);
78
 
  mpfr_set_str_raw (y, "-0.11111001001010010011010110101E-1");
79
 
  if (mpfr_cmp (z, y))
80
 
    {
81
 
      fprintf (stderr, "Error in mpfr_sub (6)\n");
82
 
      printf ("expected "); mpfr_print_binary (y); putchar ('\n');
83
 
      printf ("got      "); mpfr_print_binary (z); putchar ('\n');
84
 
      exit (1);
85
 
    }
86
 
 
87
 
  mpfr_set_prec (x, 32);
88
 
  mpfr_set_prec (y, 32);
89
 
  mpfr_set_str_raw (x, "0.10110111101001110100100101111000E0");
90
 
  mpfr_set_str_raw (y, "0.10001100100101000100110111000100E0");
91
 
  if ((inexact = mpfr_sub (x, x, y, GMP_RNDN)))
92
 
    {
93
 
      fprintf (stderr, "Wrong inexact flag (2): got %d instead of 0\n",
94
 
               inexact);
95
 
      exit (1);
96
 
    }
97
 
 
98
 
  mpfr_set_prec (x, 32);
99
 
  mpfr_set_prec (y, 32);
100
 
  mpfr_set_str_raw (x, "0.11111000110111011000100111011010E0");
101
 
  mpfr_set_str_raw (y, "0.10011111101111000100001000000000E-3");
102
 
  if ((inexact = mpfr_sub (x, x, y, GMP_RNDN)))
103
 
    {
104
 
      fprintf (stderr, "Wrong inexact flag (1): got %d instead of 0\n",
105
 
               inexact);
106
 
      exit (1);
107
 
    }
108
 
 
109
 
  mpfr_set_prec (x, 33);
110
 
  mpfr_set_prec (y, 33);
111
 
  mpfr_set_ui (x, 1, GMP_RNDN);
112
 
  mpfr_set_str_raw (y, "-0.1E-32");
113
 
  mpfr_add (x, x, y, GMP_RNDN);
114
 
  mpfr_set_str_raw (y, "0.111111111111111111111111111111111E0");
115
 
  if (mpfr_cmp (x, y))
116
 
    {
117
 
      fprintf (stderr, "Error in mpfr_sub (1 - 1E-33) with prec=33\n");
118
 
      printf ("Expected "); mpfr_print_binary (y); putchar ('\n');
119
 
      printf ("got      "); mpfr_print_binary (x); putchar ('\n');
120
 
      exit (1);
121
 
    }
122
 
 
123
 
  mpfr_set_prec (x, 32);
124
 
  mpfr_set_prec (y, 33);
125
 
  mpfr_set_ui (x, 1, GMP_RNDN);
126
 
  mpfr_set_str_raw (y, "-0.1E-32");
127
 
  mpfr_add (x, x, y, GMP_RNDN);
128
 
  if (mpfr_cmp_ui (x, 1))
129
 
    {
130
 
      fprintf (stderr, "Error in mpfr_sub (1 - 1E-33) with prec=32\n");
131
 
      printf ("Expected 1.0, got "); mpfr_print_binary (x); putchar ('\n');
132
 
      exit (1);
133
 
    }
134
 
 
135
 
  mpfr_set_prec (x, 65);
136
 
  mpfr_set_prec (y, 65);
137
 
  mpfr_set_prec (z, 64);
138
 
  mpfr_set_str_raw (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
139
 
  mpfr_set_str_raw (y, "0.1110111011110001110111011111111111101000011001011100101100101100");
140
 
  mpfr_sub (z, x, y, GMP_RNDZ);
141
 
  if (mpfr_cmp_ui (z, 1))
142
 
    {
143
 
      fprintf (stderr, "Error in mpfr_sub (1)\n");
144
 
      exit (1);
145
 
    }
146
 
  mpfr_sub (z, x, y, GMP_RNDU);
147
 
  mpfr_set_str_raw (x, "1.000000000000000000000000000000000000000000000000000000000000001");
148
 
  if (mpfr_cmp (z, x))
149
 
    {
150
 
      fprintf (stderr, "Error in mpfr_sub (2)\n");
151
 
      exit (1);
152
 
    }
153
 
  mpfr_set_str_raw (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
154
 
  mpfr_sub (z, x, y, GMP_RNDN);
155
 
  if (mpfr_cmp_ui (z, 1))
156
 
    {
157
 
      fprintf (stderr, "Error in mpfr_sub (3)\n");
158
 
      exit (1);
159
 
    }
160
 
  mpfr_set_prec (x, 66);
161
 
  mpfr_set_str_raw (x, "1.11101110111100011101110111111111111010000110010111001011001010111");
162
 
  mpfr_sub (z, x, y, GMP_RNDN);
163
 
  if (mpfr_cmp_ui (z, 1))
164
 
    {
165
 
      fprintf (stderr, "Error in mpfr_sub (4)\n");
166
 
      exit (1);
167
 
    }
168
 
  
169
 
  /* check in-place operations */
170
 
  mpfr_set_d (x, 1.0, GMP_RNDN);
171
 
  mpfr_sub (x, x, x, GMP_RNDN);
172
 
  if (mpfr_get_d1 (x) != 0.0)
173
 
    {
174
 
      fprintf (stderr, "Error for mpfr_add (x, x, x, GMP_RNDN) with x=1.0\n");
175
 
      exit (1);
176
 
    }
177
 
 
178
 
  mpfr_set_prec (x, 53);
179
 
  mpfr_set_prec (y, 53);
180
 
  mpfr_set_prec (z, 53);
181
 
  mpfr_set_d (x, 1.229318102e+09, GMP_RNDN);
182
 
  mpfr_set_d (y, 2.32221184180698677665e+05, GMP_RNDN);
183
 
  mpfr_sub (z, x, y, GMP_RNDN);
184
 
  res = 1229085880.815819263458251953125;
185
 
  got = mpfr_get_d1 (z);
186
 
  if (got != res)
187
 
    {
188
 
      fprintf (stderr, "Error in mpfr_sub (1.22e9 - 2.32e5)\n");
189
 
      printf ("expected %1.20e, got %1.20e\n", res, got);
190
 
      exit (1);
191
 
    }
192
 
 
193
 
  mpfr_set_prec (x, 112);
194
 
  mpfr_set_prec (y, 98);
195
 
  mpfr_set_prec (z, 54);
196
 
  mpfr_set_str_raw (x, "0.11111100100000000011000011100000101101010001000111E-401");
197
 
  mpfr_set_str_raw (y, "0.10110000100100000101101100011111111011101000111000101E-464");
198
 
  mpfr_sub (z, x, y, GMP_RNDN);
199
 
  if (mpfr_cmp (z, x)) {
200
 
    fprintf (stderr, "mpfr_sub(z, x, y) failed for prec(x)=112, prec(y)=98\n");
201
 
    printf ("expected "); mpfr_print_binary (x); putchar('\n');
202
 
    printf ("got      "); mpfr_print_binary (z); putchar('\n');
203
 
    exit (1);
204
 
  }
205
 
 
206
 
  mpfr_set_prec (x, 33);
207
 
  mpfr_set_ui (x, 1, GMP_RNDN);
208
 
  mpfr_div_2exp (x, x, 32, GMP_RNDN);
209
 
  mpfr_sub_ui (x, x, 1, GMP_RNDN);
210
 
 
211
 
  mpfr_set_prec (x, 5);
212
 
  mpfr_set_prec (y, 5);
213
 
  mpfr_set_str_raw (x, "1e-12");
214
 
  mpfr_set_ui (y, 1, GMP_RNDN);
215
 
  mpfr_sub (x, y, x, GMP_RNDD);
216
 
  mpfr_set_str_raw (y, "0.11111");
217
 
  if (mpfr_cmp (x, y))
218
 
    {
219
 
      fprintf (stderr, "Error in mpfr_sub (x, y, x, GMP_RNDD) for x=2^(-12), y=1\n");
220
 
      exit (1);
221
 
    }
222
 
 
223
 
  mpfr_set_prec (x, 24);
224
 
  mpfr_set_prec (y, 24);
225
 
  mpfr_set_str_raw (x, "-0.100010000000000000000000E19");
226
 
  mpfr_set_str_raw (y, "0.100000000000000000000100E15");
227
 
  mpfr_add (x, x, y, GMP_RNDD);
228
 
  mpfr_set_str_raw (y, "-0.1E19");
229
 
  if (mpfr_cmp (x, y))
230
 
    {
231
 
      fprintf (stderr, "Error in mpfr_add (2)\n");
232
 
      exit (1);
233
 
    }
234
 
 
235
 
  mpfr_set_prec (x, 2);
236
 
  mpfr_set_prec (y, 10);
237
 
  mpfr_set_prec (z, 10);
238
 
  mpfr_set_ui (y, 0, GMP_RNDN);
239
 
  mpfr_set_str_raw (z, "0.100000000000000000000100E15");
240
 
  if (mpfr_sub (x, y, z, GMP_RNDN) <= 0)
241
 
    {
242
 
      fprintf (stderr, "Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n");
243
 
      exit (1);
244
 
    }
245
 
  if (mpfr_sub (x, z, y, GMP_RNDN) >= 0)
246
 
    {
247
 
      fprintf (stderr, "Wrong inexact flag in x=mpfr_sub(z,0) for prec(z)>prec(x)\n");
248
 
      exit (1);
249
 
    }
250
 
 
251
 
  mpfr_clear (x);
252
 
  mpfr_clear (y);
253
 
  mpfr_clear (z);
254
 
}
255
 
 
256
 
void
257
 
bug_ddefour(void)
258
 
{
259
 
    mpfr_t ex, ex1, ex2, ex3, tot, tot1;
260
 
 
261
 
    mpfr_init2(ex, 53);
262
 
    mpfr_init2(ex1, 53);
263
 
    mpfr_init2(ex2, 53);
264
 
    mpfr_init2(ex3, 53);
265
 
    mpfr_init2(tot, 150);
266
 
    mpfr_init2(tot1, 150);
267
 
 
268
 
    mpfr_set_ui( ex, 1, GMP_RNDN);
269
 
    mpfr_mul_2exp( ex, ex, 906, GMP_RNDN);
270
 
    mpfr_log( tot, ex, GMP_RNDN);
271
 
    mpfr_set( ex1, tot, GMP_RNDN); /* ex1 = high(tot) */
272
 
    mpfr_sub( ex2, tot, ex1, GMP_RNDN); /* ex2 = high(tot - ex1) */
273
 
    mpfr_sub( tot1, tot, ex1, GMP_RNDN); /* tot1 = tot - ex1 */
274
 
    mpfr_set( ex3, tot1, GMP_RNDN); /* ex3 = high(tot - ex1) */
275
 
 
276
 
    if (mpfr_cmp(ex2, ex3)) 
277
 
      {
278
 
        fprintf (stderr, "Error in ddefour test.\n");
279
 
        printf ("ex2="); mpfr_print_binary (ex2); putchar ('\n');
280
 
        printf ("ex3="); mpfr_print_binary (ex3); putchar ('\n');
281
 
        exit (1);
282
 
      }
283
 
 
284
 
    mpfr_clear (ex);
285
 
    mpfr_clear (ex1);
286
 
    mpfr_clear (ex2);
287
 
    mpfr_clear (ex3);
288
 
    mpfr_clear (tot);
289
 
    mpfr_clear (tot1);
290
 
}
291
 
 
292
 
/* if u = o(x-y), v = o(u-x), w = o(v+y), then x-y = u-w */
293
 
void
294
 
check_two_sum (mp_prec_t p)
295
 
{
296
 
  mpfr_t x, y, u, v, w;
297
 
  mp_rnd_t rnd;
298
 
  int inexact;
299
 
  
300
 
  mpfr_init2 (x, p);
301
 
  mpfr_init2 (y, p);
302
 
  mpfr_init2 (u, p);
303
 
  mpfr_init2 (v, p);
304
 
  mpfr_init2 (w, p);
305
 
  mpfr_random (x);
306
 
  mpfr_random (y);
307
 
  if (mpfr_cmp_abs (x, y) < 0)
308
 
    mpfr_swap (x, y);
309
 
  rnd = LONG_RAND() % 4;
310
 
  rnd = GMP_RNDN;
311
 
  inexact = mpfr_sub (u, x, y, GMP_RNDN);
312
 
  mpfr_sub (v, u, x, GMP_RNDN);
313
 
  mpfr_add (w, v, y, GMP_RNDN);
314
 
  /* as u = (x-y) - w, we should have inexact and w of opposite signs */
315
 
  if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
316
 
      ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) ||
317
 
      ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0)))
318
 
    {
319
 
      fprintf (stderr, "Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p,
320
 
               mpfr_print_rnd_mode (rnd));
321
 
      printf ("x="); mpfr_print_binary(x); putchar('\n');
322
 
      printf ("y="); mpfr_print_binary(y); putchar('\n');
323
 
      printf ("u="); mpfr_print_binary(u); putchar('\n');
324
 
      printf ("v="); mpfr_print_binary(v); putchar('\n');
325
 
      printf ("w="); mpfr_print_binary(w); putchar('\n');
326
 
      printf ("inexact = %d\n", inexact);
327
 
      exit (1);
328
 
    }
329
 
  mpfr_clear (x);
330
 
  mpfr_clear (y);
331
 
  mpfr_clear (u);
332
 
  mpfr_clear (v);
333
 
  mpfr_clear (w);
334
 
}
335
 
 
336
 
#define MAX_PREC 100
337
 
 
338
 
void
339
 
check_inexact (void)
340
 
{
341
 
  mpfr_t x, y, z, u;
342
 
  mp_prec_t px, py, pu, pz;
343
 
  int inexact, cmp;
344
 
  mp_rnd_t rnd;
345
 
  
346
 
  mpfr_init (x);
347
 
  mpfr_init (y);
348
 
  mpfr_init (z);
349
 
  mpfr_init (u);
350
 
 
351
 
  for (px=2; px<MAX_PREC; px++)
352
 
    {
353
 
      mpfr_set_prec (x, px);
354
 
      mpfr_random (x);
355
 
      for (pu=2; pu<MAX_PREC; pu++)
356
 
        {
357
 
          mpfr_set_prec (u, pu);
358
 
          mpfr_random (u);
359
 
          for (py=2; py<MAX_PREC; py++)
360
 
            {
361
 
              mpfr_set_prec (y, py);
362
 
              pz =  (mpfr_cmp_abs (x, u) >= 0) ? MPFR_EXP(x)-MPFR_EXP(u)
363
 
                : MPFR_EXP(u)-MPFR_EXP(x);
364
 
              pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u));
365
 
              mpfr_set_prec (z, pz);
366
 
              rnd = LONG_RAND () % 4;
367
 
              if (mpfr_sub (z, x, u, rnd))
368
 
                {
369
 
                  fprintf (stderr, "z <- x - u should be exact\n");
370
 
                  exit (1);
371
 
                }
372
 
              for (rnd=0; rnd<4; rnd++)
373
 
                {
374
 
                  inexact = mpfr_sub (y, x, u, rnd);
375
 
                  cmp = mpfr_cmp (y, z);
376
 
                  if (((inexact == 0) && (cmp != 0)) ||
377
 
                      ((inexact > 0) && (cmp <= 0)) ||
378
 
                      ((inexact < 0) && (cmp >= 0)))
379
 
                    {
380
 
                      fprintf (stderr, "Wrong inexact flag for rnd=%s\n",
381
 
                           mpfr_print_rnd_mode(rnd));
382
 
                      printf ("expected %d, got %d\n", cmp, inexact);
383
 
                      printf ("x="); mpfr_print_binary (x); putchar ('\n');
384
 
                      printf ("u="); mpfr_print_binary (u); putchar ('\n');
385
 
                      printf ("y=  "); mpfr_print_binary (y); putchar ('\n');
386
 
                      printf ("x-u="); mpfr_print_binary (z); putchar ('\n');
387
 
                      exit (1);
388
 
                    }
389
 
                }
390
 
            }
391
 
        }
392
 
    }
393
 
 
394
 
  mpfr_clear (x);
395
 
  mpfr_clear (y);
396
 
  mpfr_clear (z);
397
 
  mpfr_clear (u);
398
 
}
399
 
 
400
 
int
401
 
main(void)
402
 
{
403
 
  mp_prec_t p;
404
 
  unsigned i;
405
 
 
406
 
  check_diverse ();
407
 
  check_inexact ();
408
 
  bug_ddefour ();
409
 
 
410
 
  for (p=2; p<200; p++)
411
 
    for (i=0; i<200; i++)
412
 
      check_two_sum (p);
413
 
 
414
 
  return 0;
415
 
}