1
/* Test file for mpfr_add and mpfr_sub.
3
Copyright 1999, 2000, 2001, 2002 Free Software Foundation.
5
This file is part of the MPFR Library.
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.
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.
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. */
29
#include "mpfr-impl.h"
30
#include "mpfr-test.h"
32
void checknan _PROTO((double, double, mp_rnd_t, unsigned int, unsigned int, unsigned int));
33
void check3 _PROTO((double, double, mp_rnd_t));
34
void check4 _PROTO((double, double, mp_rnd_t));
35
void check5 _PROTO((double, mp_rnd_t));
36
void check2 _PROTO((double, int, double, int, int, int));
37
void check2a _PROTO((double, int, double, int, int, int, char *));
38
void check64 _PROTO((void));
39
void check_same _PROTO((void));
40
void check_case_1b _PROTO((void));
41
void check_case_2 _PROTO((void));
42
void check_inexact _PROTO((void));
45
/* Parameter "z1" of check() used to be last in the argument list, but that
46
tickled a bug in 32-bit sparc gcc 2.95.2. A "double" in that position is
47
passed on the stack at an address which is 4mod8, but the generated code
48
didn't take into account that alignment, resulting in bus errors. The
49
easiest workaround is to move it to the start of the arg list (where it's
50
passed in registers), this macro does that. FIXME: Change the actual
51
calls to check(), rather than using a macro. */
53
#define check(x,y,rnd_mode,px,py,pz,z1) _check(x,y,z1,rnd_mode,px,py,pz)
55
/* checks that x+y gives the same results in double
56
and with mpfr with 53 bits of precision */
58
_check (double x, double y, double z1, mp_rnd_t rnd_mode, unsigned int px,
59
unsigned int py, unsigned int pz)
61
double z2; mpfr_t xx,yy,zz; int cert=0;
66
mpfr_set_d(xx, x, rnd_mode);
67
mpfr_set_d(yy, y, rnd_mode);
68
mpfr_add(zz, xx, yy, rnd_mode);
69
#ifdef MPFR_HAVE_FESETROUND
70
mpfr_set_machine_rnd_mode(rnd_mode);
71
if (px==53 && py==53 && pz==53) cert=1;
73
if (z1==0.0) z1=x+y; else cert=1;
74
z2 = mpfr_get_d1 (zz);
75
mpfr_set_d (yy, z2, GMP_RNDN);
76
if (!mpfr_cmp (zz, yy) && cert && z1!=z2 && !(isnan(z1) && isnan(z2))) {
77
printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
78
printf("mpfr_add failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
79
x, y, mpfr_print_rnd_mode(rnd_mode));
82
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
86
checknan (double x, double y, mp_rnd_t rnd_mode, unsigned int px,
87
unsigned int py, unsigned int pz)
89
double z2; mpfr_t xx, yy, zz;
94
mpfr_set_d(xx, x, rnd_mode);
95
mpfr_set_d(yy, y, rnd_mode);
96
mpfr_add(zz, xx, yy, rnd_mode);
97
#ifdef MPFR_HAVE_FESETROUND
98
mpfr_set_machine_rnd_mode(rnd_mode);
100
if (MPFR_IS_NAN(zz) == 0) { printf("Error, not an MPFR_NAN for xx = %1.20e, y = %1.20e\n", x, y); exit(1); }
101
z2 = mpfr_get_d1 (zz);
102
if (!isnan(z2)) { printf("Error, not a NaN after conversion, xx = %1.20e yy = %1.20e, got %1.20e\n", x, y, z2); exit(1); }
104
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
107
#ifdef MPFR_HAVE_FESETROUND
108
/* idem than check for mpfr_add(x, x, y) */
110
check3 (double x, double y, mp_rnd_t rnd_mode)
112
double z1,z2; mpfr_t xx,yy; int neg;
114
neg = LONG_RAND() % 2;
117
mpfr_set_d(xx, x, rnd_mode);
118
mpfr_set_d(yy, y, rnd_mode);
119
if (neg) mpfr_sub(xx, xx, yy, rnd_mode);
120
else mpfr_add(xx, xx, yy, rnd_mode);
121
mpfr_set_machine_rnd_mode(rnd_mode);
122
z1 = (neg) ? x-y : x+y;
123
z2 = mpfr_get_d1 (xx);
124
mpfr_set_d (yy, z2, GMP_RNDN);
125
if (!mpfr_cmp (xx, yy) && z1!=z2 && !(isnan(z1) && isnan(z2))) {
126
printf("expected result is %1.20e, got %1.20e\n",z1,z2);
127
printf("mpfr_%s(x,x,y) failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",
128
(neg) ? "sub" : "add",x,y,rnd_mode);
131
mpfr_clear(xx); mpfr_clear(yy);
134
/* idem than check for mpfr_add(x, y, x) */
136
check4 (double x, double y, mp_rnd_t rnd_mode)
142
neg = LONG_RAND() % 2;
145
mpfr_set_d(xx, x, rnd_mode);
146
mpfr_set_d(yy, y, rnd_mode);
147
if (neg) mpfr_sub(xx, yy, xx, rnd_mode);
148
else mpfr_add(xx, yy, xx, rnd_mode);
149
mpfr_set_machine_rnd_mode(rnd_mode);
150
z1 = (neg) ? y-x : x+y;
151
z2 = mpfr_get_d1 (xx);
152
mpfr_set_d (yy, z2, GMP_RNDN);
153
/* check that xx is representable as a double and no overflow occurred */
154
if ((mpfr_cmp (xx, yy) == 0) && (z1 != z2)) {
155
printf("expected result is %1.20e, got %1.20e\n", z1, z2);
156
printf("mpfr_%s(x,y,x) failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
157
(neg) ? "sub" : "add", x, y, mpfr_print_rnd_mode(rnd_mode));
160
mpfr_clear(xx); mpfr_clear(yy);
163
/* idem than check for mpfr_add(x, x, x) */
165
check5 (double x, mp_rnd_t rnd_mode)
167
double z1,z2; mpfr_t xx, yy; int neg;
171
neg = LONG_RAND() % 2;
172
mpfr_set_d(xx, x, rnd_mode);
173
if (neg) mpfr_sub(xx, xx, xx, rnd_mode);
174
else mpfr_add(xx, xx, xx, rnd_mode);
175
mpfr_set_machine_rnd_mode(rnd_mode);
176
z1 = (neg) ? x-x : x+x;
177
z2 = mpfr_get_d1 (xx);
178
mpfr_set_d (yy, z2, GMP_RNDN);
179
/* check NaNs first since mpfr_cmp does not like them */
180
if (!(isnan(z1) && isnan(z2)) && !mpfr_cmp (xx, yy) && z1!=z2)
182
printf ("expected result is %1.20e, got %1.20e\n",z1,z2);
183
printf ("mpfr_%s(x,x,x) failed for x=%1.20e with rnd_mode=%s\n",
184
(neg) ? "sub" : "add", x, mpfr_print_rnd_mode (rnd_mode));
192
check2 (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode)
194
mpfr_t xx, yy, zz; double z,z2; int u;
196
mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
197
mpfr_set_d(xx, x, rnd_mode);
198
mpfr_set_d(yy, y, rnd_mode);
199
mpfr_add(zz, xx, yy, rnd_mode);
200
mpfr_set_machine_rnd_mode(rnd_mode);
201
z = x+y; z2=mpfr_get_d1 (zz); u=ulp(z,z2);
202
/* one ulp difference is possible due to composed rounding */
203
if (px>=53 && py>=53 && pz>=53 && ABS(u)>1) {
204
printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
205
x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
206
printf("got %1.20e\n",z2);
207
printf("result should be %1.20e (diff=%d ulp)\n",z,u);
208
mpfr_set_d(zz, z, rnd_mode);
209
printf("i.e."); mpfr_print_binary(zz); putchar('\n');
211
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
216
check2a (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode,
221
mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
222
mpfr_set_d(xx, x, rnd_mode);
223
mpfr_set_d(yy, y, rnd_mode);
224
mpfr_add(zz, xx, yy, rnd_mode);
225
mpfr_set_prec(xx, pz);
226
mpfr_set_str(xx, res, 16, GMP_RNDN);
227
if (mpfr_cmp(xx, zz)) {
228
printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
229
x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
230
printf("got "); mpfr_print_binary(zz); putchar('\n');
231
printf("instead of "); mpfr_print_binary(xx); putchar('\n');
234
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
246
mpfr_set_prec (x, 29);
247
mpfr_set_str_raw (x, "1.1101001000101111011010010110e-3");
248
mpfr_set_prec (t, 58);
249
mpfr_set_str_raw (t, "0.11100010011111001001100110010111110110011000000100101E-1");
250
mpfr_set_prec (u, 29);
251
mpfr_add (u, x, t, GMP_RNDD);
252
mpfr_set_str_raw (t, "1.0101011100001000011100111110e-1");
255
fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=29, prec(t)=58\n");
256
printf ("expected "); mpfr_out_str (stdout, 2, 29, t, GMP_RNDN);
258
printf ("got "); mpfr_out_str (stdout, 2, 29, u, GMP_RNDN);
263
mpfr_set_prec (x, 4);
264
mpfr_set_str_raw (x, "-1.0E-2");
265
mpfr_set_prec (t, 2);
266
mpfr_set_str_raw (t, "-1.1e-2");
267
mpfr_set_prec (u, 2);
268
mpfr_add (u, x, t, GMP_RNDN);
269
if (MPFR_MANT(u)[0] << 2)
271
fprintf (stderr, "result not normalized for prec=2\n");
272
mpfr_print_binary (u); putchar ('\n');
275
mpfr_set_str_raw (t, "-1.0e-1");
278
fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=4, prec(t)=2\n");
279
printf ("expected -1.0e-1\n");
280
printf ("got "); mpfr_out_str (stdout, 2, 4, u, GMP_RNDN);
285
mpfr_set_prec (x, 8);
286
mpfr_set_str_raw (x, "-0.10011010"); /* -77/128 */
287
mpfr_set_prec (t, 4);
288
mpfr_set_str_raw (t, "-1.110e-5"); /* -7/128 */
289
mpfr_set_prec (u, 4);
290
mpfr_add (u, x, t, GMP_RNDN); /* should give -5/8 */
291
mpfr_set_str_raw (t, "-1.010e-1");
292
if (mpfr_cmp (u, t)) {
293
fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=8, prec(t)=4\n");
294
printf ("expected -1.010e-1\n");
295
printf ("got "); mpfr_out_str (stdout, 2, 4, u, GMP_RNDN);
300
mpfr_set_prec (x, 112); mpfr_set_prec (t, 98); mpfr_set_prec (u, 54);
301
mpfr_set_str_raw (x, "-0.11111100100000000011000011100000101101010001000111E-401");
302
mpfr_set_str_raw (t, "0.10110000100100000101101100011111111011101000111000101E-464");
303
mpfr_add (u, x, t, GMP_RNDN);
304
if (mpfr_cmp (u, x)) {
305
fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=112, prec(t)=98\n");
309
mpfr_set_prec (x, 92); mpfr_set_prec (t, 86); mpfr_set_prec (u, 53);
310
mpfr_set_d (x, -5.03525136761487735093e-74, GMP_RNDN);
311
mpfr_set_d (t, 8.51539046314262304109e-91, GMP_RNDN);
312
mpfr_add (u, x, t, GMP_RNDN);
313
if (mpfr_get_d1 (u) != -5.0352513676148773509283672e-74) {
314
fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=92, prec(t)=86\n");
318
mpfr_set_prec(x, 53); mpfr_set_prec(t, 76); mpfr_set_prec(u, 76);
319
mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
320
mpfr_set_str_raw(t, "-0.1011000101110010000101111111011111010001110011110111100110101011110010011111");
321
mpfr_sub(u, x, t, GMP_RNDU);
322
mpfr_set_str_raw(t, "0.1011000101110010000101111111011100111111101010011011110110101011101000000100");
324
printf("expect "); mpfr_print_binary(t); putchar('\n');
325
fprintf (stderr, "mpfr_add failed for precisions 53-76\n"); exit(1);
327
mpfr_set_prec(x, 53); mpfr_set_prec(t, 108); mpfr_set_prec(u, 108);
328
mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
329
mpfr_set_str_raw(t, "-0.101100010111001000010111111101111101000111001111011110011010101111001001111000111011001110011000000000111111");
330
mpfr_sub(u, x, t, GMP_RNDU);
331
mpfr_set_str_raw(t, "0.101100010111001000010111111101110011111110101001101111011010101110100000001011000010101110011000000000111111");
333
printf("expect "); mpfr_print_binary(t); putchar('\n');
334
fprintf(stderr, "mpfr_add failed for precisions 53-108\n"); exit(1);
336
mpfr_set_prec(x, 97); mpfr_set_prec(t, 97); mpfr_set_prec(u, 97);
337
mpfr_set_str_raw(x, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010000000000000000E-39");
338
mpfr_set_ui(t, 1, GMP_RNDN);
339
mpfr_add(u, x, t, GMP_RNDN);
340
mpfr_set_str_raw(x, "0.1000000000000000000000000000000000000000111110110000100000000101100011011110100000101111100010001E1");
342
fprintf(stderr, "mpfr_add failed for precision 97\n"); exit(1);
344
mpfr_set_prec(x, 128); mpfr_set_prec(t, 128); mpfr_set_prec(u, 128);
345
mpfr_set_str_raw(x, "0.10101011111001001010111011001000101100111101000000111111111011010100001100011101010001010111111101111010100110111111100101100010E-4");
346
mpfr_set(t, x, GMP_RNDN);
347
mpfr_sub(u, x, t, GMP_RNDN);
348
mpfr_set_prec(x, 96); mpfr_set_prec(t, 96); mpfr_set_prec(u, 96);
349
mpfr_set_str_raw(x, "0.111000000001110100111100110101101001001010010011010011100111100011010100011001010011011011000010E-4");
350
mpfr_set(t, x, GMP_RNDN);
351
mpfr_sub(u, x, t, GMP_RNDN);
352
mpfr_set_prec(x, 85); mpfr_set_prec(t, 85); mpfr_set_prec(u, 85);
353
mpfr_set_str_raw(x, "0.1111101110100110110110100010101011101001100010100011110110110010010011101100101111100E-4");
354
mpfr_set_str_raw(t, "0.1111101110100110110110100010101001001000011000111000011101100101110100001110101010110E-4");
355
mpfr_sub(u, x, t, GMP_RNDU);
356
mpfr_sub(x, x, t, GMP_RNDU);
357
if (mpfr_cmp(x, u) != 0) {
358
printf("Error in mpfr_sub: u=x-t and x=x-t give different results\n");
361
if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
362
((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
363
printf("Error in mpfr_sub: result is not msb-normalized (1)\n"); exit(1);
365
mpfr_set_prec(x, 65); mpfr_set_prec(t, 65); mpfr_set_prec(u, 65);
366
mpfr_set_str_raw(x, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
367
mpfr_set_str_raw(t, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
368
mpfr_sub(u, x, t, GMP_RNDU);
369
if (mpfr_get_d1 (u) != 9.4349060620538533806e167) { /* 2^558 */
370
printf("Error (1) in mpfr_sub\n"); exit(1);
373
mpfr_set_prec(x, 64); mpfr_set_prec(t, 64); mpfr_set_prec(u, 64);
374
mpfr_set_str_raw(x, "0.1000011110101111011110111111000011101011101111101101101100000100E-220");
375
mpfr_set_str_raw(t, "0.1000011110101111011110111111000011101011101111101101010011111101E-220");
376
mpfr_add(u, x, t, GMP_RNDU);
377
if ((MPFR_MANT(u)[0] & 1) != 1) {
378
printf("error in mpfr_add with rnd_mode=GMP_RNDU\n");
379
printf("b= "); mpfr_print_binary(x); putchar('\n');
380
printf("c= "); mpfr_print_binary(t); putchar('\n');
381
printf("b+c="); mpfr_print_binary(u); putchar('\n');
385
/* bug found by Norbert Mueller, 14 Sep 2000 */
386
mpfr_set_prec(x, 56); mpfr_set_prec(t, 83); mpfr_set_prec(u, 10);
387
mpfr_set_str_raw(x, "0.10001001011011001111101100110100000101111010010111010111E-7");
388
mpfr_set_str_raw(t, "0.10001001011011001111101100110100000101111010010111010111000000000111110110110000100E-7");
389
mpfr_sub(u, x, t, GMP_RNDU);
391
/* array bound write found by Norbert Mueller, 26 Sep 2000 */
392
mpfr_set_prec(x, 109); mpfr_set_prec(t, 153); mpfr_set_prec(u, 95);
393
mpfr_set_str_raw(x,"0.1001010000101011101100111000110001111111111111111111111111111111111111111111111111111111111111100000000000000E33");
394
mpfr_set_str_raw(t,"-0.100101000010101110110011100011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011100101101000000100100001100110111E33");
395
mpfr_add(u, x, t, GMP_RNDN);
397
/* array bound writes found by Norbert Mueller, 27 Sep 2000 */
398
mpfr_set_prec(x, 106); mpfr_set_prec(t, 53); mpfr_set_prec(u, 23);
399
mpfr_set_str_raw(x, "-0.1000011110101111111001010001000100001011000000000000000000000000000000000000000000000000000000000000000000E-59");
400
mpfr_set_str_raw(t, "-0.10000111101011111110010100010001101100011100110100000E-59");
401
mpfr_sub(u, x, t, GMP_RNDN);
402
mpfr_set_prec(x, 177); mpfr_set_prec(t, 217); mpfr_set_prec(u, 160);
403
mpfr_set_str_raw(x, "-0.111010001011010000111001001010010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E35");
404
mpfr_set_str_raw(t, "0.1110100010110100001110010010100100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111011010011100001111001E35");
405
mpfr_add(u, x, t, GMP_RNDN);
406
mpfr_set_prec(x, 214); mpfr_set_prec(t, 278); mpfr_set_prec(u, 207);
407
mpfr_set_str_raw(x, "0.1000100110100110101101101101000000010000100111000001001110001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E66");
408
mpfr_set_str_raw(t, "-0.10001001101001101011011011010000000100001001110000010011100010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001111011111001001100011E66");
409
mpfr_add(u, x, t, GMP_RNDN);
410
mpfr_set_prec(x, 32); mpfr_set_prec(t, 247); mpfr_set_prec(u, 223);
411
mpfr_set_str_raw(x, "0.10000000000000000000000000000000E1");
412
mpfr_set_str_raw(t, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100000100011110000101110110011101110100110110111111011010111100100000000000000000000000000E0");
413
mpfr_sub(u, x, t, GMP_RNDN);
414
if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
415
((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
416
printf("Error in mpfr_sub: result is not msb-normalized (2)\n"); exit(1);
419
/* bug found by Nathalie Revol, 21 March 2001 */
420
mpfr_set_prec (x, 65);
421
mpfr_set_prec (t, 65);
422
mpfr_set_prec (u, 65);
423
mpfr_set_str_raw (x, "0.11100100101101001100111011111111110001101001000011101001001010010E-35");
424
mpfr_set_str_raw (t, "0.10000000000000000000000000000000000001110010010110100110011110000E1");
425
mpfr_sub (u, t, x, GMP_RNDU);
426
if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
427
((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
428
fprintf(stderr, "Error in mpfr_sub: result is not msb-normalized (3)\n");
432
/* bug found by Fabrice Rouillier, 27 Mar 2001 */
433
mpfr_set_prec (x, 107);
434
mpfr_set_prec (t, 107);
435
mpfr_set_prec (u, 107);
436
mpfr_set_str_raw (x, "0.10111001001111010010001000000010111111011011011101000001001000101000000000000000000000000000000000000000000E315");
437
mpfr_set_str_raw (t, "0.10000000000000000000000000000000000101110100100101110110000001100101011111001000011101111100100100111011000E350");
438
mpfr_sub (u, x, t, GMP_RNDU);
439
if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
440
((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
441
fprintf(stderr, "Error in mpfr_sub: result is not msb-normalized (4)\n");
445
/* checks that NaN flag is correctly reset */
446
mpfr_set_d (t, 1.0, GMP_RNDN);
447
mpfr_set_d (u, 1.0, GMP_RNDN);
449
mpfr_add (x, t, u, GMP_RNDN);
450
if (mpfr_cmp_ui (x, 2)) {
451
fprintf (stderr, "Error in mpfr_add: 1+1 gives %e\n", mpfr_get_d1 (x));
455
mpfr_clear(x); mpfr_clear(t); mpfr_clear(u);
458
/* check case when c does not overlap with a, but both b and c count
464
unsigned int prec_a, prec_b, prec_c, dif;
470
for (prec_a = 2; prec_a <= 64; prec_a++)
472
mpfr_set_prec (a, prec_a);
473
for (prec_b = prec_a + 2; prec_b <= 64; prec_b++)
475
dif = prec_b - prec_a;
476
mpfr_set_prec (b, prec_b);
477
/* b = 1 - 2^(-prec_a) + 2^(-prec_b) */
478
mpfr_set_ui (b, 1, GMP_RNDN);
479
mpfr_div_2exp (b, b, dif, GMP_RNDN);
480
mpfr_sub_ui (b, b, 1, GMP_RNDN);
481
mpfr_div_2exp (b, b, prec_a, GMP_RNDN);
482
mpfr_add_ui (b, b, 1, GMP_RNDN);
483
for (prec_c = dif; prec_c <= 64; prec_c++)
485
/* c = 2^(-prec_a) - 2^(-prec_b) */
486
mpfr_set_prec (c, prec_c);
487
mpfr_set_si (c, -1, GMP_RNDN);
488
mpfr_div_2exp (c, c, dif, GMP_RNDN);
489
mpfr_add_ui (c, c, 1, GMP_RNDN);
490
mpfr_div_2exp (c, c, prec_a, GMP_RNDN);
491
mpfr_add (a, b, c, GMP_RNDN);
492
if (mpfr_cmp_ui (a, 1) != 0)
494
fprintf (stderr, "case (1b) failed for prec_a=%u, prec_b=%u, prec_c=%u\n", prec_a, prec_b, prec_c);
495
printf("b="); mpfr_print_binary(b); putchar('\n');
496
printf("c="); mpfr_print_binary(c); putchar('\n');
497
printf("a="); mpfr_print_binary(a); putchar('\n');
509
/* check case when c overlaps with a */
520
mpfr_set_str_raw(a, "1E110"); /* a = 2^110 */
521
mpfr_set_str_raw(b, "1E900"); /* b = 2^900 */
522
mpfr_set_str_raw(c, "1E500"); /* c = 2^500 */
523
mpfr_add(c, c, a, GMP_RNDZ); /* c = 2^500 + 2^110 */
524
mpfr_sub(d, b, c, GMP_RNDZ); /* d = 2^900 - 2^500 - 2^110 */
525
mpfr_add(b, b, c, GMP_RNDZ); /* b = 2^900 + 2^500 + 2^110 */
526
mpfr_add(a, b, d, GMP_RNDZ); /* a = 2^901 */
527
if (mpfr_cmp_ui_2exp (a, 1, 901))
529
fprintf (stderr, "b + d fails for b=2^900+2^500+2^110, d=2^900-2^500-2^110\n");
530
fprintf (stderr, "expected 1.0e901, got ");
531
mpfr_out_str (stderr, 2, 0, a, GMP_RNDN);
532
fprintf (stderr, "\n");
542
/* checks when source and destination are equal */
548
mpfr_init(x); mpfr_set_d(x, 1.0, GMP_RNDZ);
549
mpfr_add(x, x, x, GMP_RNDZ);
550
if (mpfr_get_d1 (x) != 2.0) {
551
printf("Error when all 3 operands are equal\n"); exit(1);
556
#define check53(x, y, r, z) check(x, y, r, 53, 53, 53, z)
557
#define check53nan(x, y, r) checknan(x, y, r, 53, 53, 53);
565
mp_prec_t px, py, pu, pz;
574
mpfr_set_prec (x, 2);
575
mpfr_set_str_raw (x, "0.1E-4");
576
mpfr_set_prec (u, 33);
577
mpfr_set_str_raw (u, "0.101110100101101100000000111100000E-1");
578
mpfr_set_prec (y, 31);
579
if ((inexact = mpfr_add (y, x, u, GMP_RNDN)))
581
fprintf (stderr, "Wrong inexact flag (2): expected 0, got %d\n", inexact);
585
mpfr_set_prec (x, 2);
586
mpfr_set_str_raw (x, "0.1E-4");
587
mpfr_set_prec (u, 33);
588
mpfr_set_str_raw (u, "0.101110100101101100000000111100000E-1");
589
mpfr_set_prec (y, 28);
590
if ((inexact = mpfr_add (y, x, u, GMP_RNDN)))
592
fprintf (stderr, "Wrong inexact flag (1): expected 0, got %d\n", inexact);
596
for (px=2; px<MAX_PREC; px++)
598
mpfr_set_prec (x, px);
600
for (pu=2; pu<MAX_PREC; pu++)
602
mpfr_set_prec (u, pu);
604
for (py=2; py<MAX_PREC; py++)
606
mpfr_set_prec (y, py);
607
pz = (mpfr_cmp_abs (x, u) >= 0) ? MPFR_EXP(x)-MPFR_EXP(u)
608
: MPFR_EXP(u)-MPFR_EXP(x);
609
/* x + u is exactly representable with precision
610
abs(EXP(x)-EXP(u)) + max(prec(x), prec(u)) + 1 */
611
pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u)) + 1;
612
mpfr_set_prec (z, pz);
613
rnd = LONG_RAND () % 4;
614
if (mpfr_add (z, x, u, rnd))
616
fprintf (stderr, "z <- x + u should be exact\n");
617
printf ("x="); mpfr_print_binary (x); putchar ('\n');
618
printf ("u="); mpfr_print_binary (u); putchar ('\n');
619
printf ("z="); mpfr_print_binary (z); putchar ('\n');
622
for (rnd=0; rnd<4; rnd++)
624
inexact = mpfr_add (y, x, u, rnd);
625
cmp = mpfr_cmp (y, z);
626
if (((inexact == 0) && (cmp != 0)) ||
627
((inexact > 0) && (cmp <= 0)) ||
628
((inexact < 0) && (cmp >= 0)))
630
fprintf (stderr, "Wrong inexact flag for rnd=%s\n",
631
mpfr_print_rnd_mode(rnd));
632
printf ("expected %d, got %d\n", cmp, inexact);
633
printf ("x="); mpfr_print_binary (x); putchar ('\n');
634
printf ("u="); mpfr_print_binary (u); putchar ('\n');
635
printf ("y= "); mpfr_print_binary (y); putchar ('\n');
636
printf ("x+u="); mpfr_print_binary (z); putchar ('\n');
651
main (int argc, char *argv[])
653
#ifdef MPFR_HAVE_FESETROUND
666
check(293607738.0, 1.9967571564050541e-5, GMP_RNDU, 64, 53, 53,
667
2.9360773800002003e8);
668
check(880524.0, -2.0769715792901673e-5, GMP_RNDN, 64, 53, 53,
669
8.8052399997923023e5);
670
check(1196426492.0, -1.4218093058435347e-3, GMP_RNDN, 64, 53, 53,
671
1.1964264919985781e9);
672
check(982013018.0, -8.941829477291838e-7, GMP_RNDN, 64, 53, 53,
673
9.8201301799999905e8);
674
check(1092583421.0, 1.0880649218158844e9, GMP_RNDN, 64, 53, 53,
675
2.1806483428158846e9);
676
check(1.8476886419022969e-6, 961494401.0, GMP_RNDN, 53, 64, 53,
677
9.6149440100000179e8);
678
check(-2.3222118418069868e5, 1229318102.0, GMP_RNDN, 53, 64, 53,
679
1.2290858808158193e9);
680
check(-3.0399171300395734e-6, 874924868.0, GMP_RNDN, 53, 64, 53,
681
8.749248679999969e8);
682
check(9.064246624706179e1, 663787413.0, GMP_RNDN, 53, 64, 53,
683
6.6378750364246619e8);
684
check(-1.0954322421551264e2, 281806592.0, GMP_RNDD, 53, 64, 53,
685
2.8180648245677572e8);
686
check(5.9836930386056659e-8, 1016217213.0, GMP_RNDN, 53, 64, 53,
687
1.0162172130000001e9);
688
check(-1.2772161928500301e-7, 1237734238.0, GMP_RNDN, 53, 64, 53,
689
1.2377342379999998e9);
690
check(-4.567291988483277e8, 1262857194.0, GMP_RNDN, 53, 64, 53,
691
8.0612799515167236e8);
692
check(4.7719471752925262e7, 196089880.0, GMP_RNDN, 53, 53, 53,
693
2.4380935175292528e8);
694
check(4.7719471752925262e7, 196089880.0, GMP_RNDN, 53, 64, 53,
695
2.4380935175292528e8);
696
check(-1.716113812768534e-140, 1271212614.0, GMP_RNDZ, 53, 64, 53,
697
1.2712126139999998e9);
698
check(-1.2927455200185474e-50, 1675676122.0, GMP_RNDD, 53, 64, 53,
699
1.6756761219999998e9);
700
check53(1.22191250737771397120e+20, 948002822.0, GMP_RNDN,
701
122191250738719408128.0);
702
check53(9966027674114492.0, 1780341389094537.0, GMP_RNDN,
703
11746369063209028.0);
704
check53(2.99280481918991653800e+272, 5.34637717585790933424e+271, GMP_RNDN,
705
3.5274425367757071711e272);
707
check53(6.14384195492641560499e-02, -6.14384195401037683237e-02, GMP_RNDU,
708
9.1603877261370314499e-12);
709
check53(1.16809465359248765399e+196, 7.92883212101990665259e+196, GMP_RNDU,
710
9.0969267746123943065e196);
711
check53(3.14553393112021279444e-67, 3.14553401015952024126e-67, GMP_RNDU,
712
6.2910679412797336946e-67);
714
SEED_RAND (time(NULL));
715
check53(5.43885304644369509058e+185,-1.87427265794105342763e-57,GMP_RNDN,
716
5.4388530464436950905e185);
717
check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDZ,
718
5.4388530464436944867e185);
719
check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDU,
720
5.4388530464436950905e185);
721
check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDD,
722
5.4388530464436944867e185);
723
check2a(6.85523243386777784171e+107,187,-2.78148588123699111146e+48,87,178,
724
GMP_RNDD, "4.ab980a5cb9407ffffffffffffffffffffffffffffffe@89");
725
check2a(-1.21510626304662318398e+145,70,1.21367733647758957118e+145,65,61,
726
GMP_RNDD, "-1.2bfad031d94@118");
727
check2a(2.73028857032080744543e+155,83,-1.16446121423113355603e+163,59,125,
728
GMP_RNDZ, "-3.3c42dee09703d0639a6@135");
729
check2a(-4.38589520019641698848e+78,155,-1.09923643769309483415e+72,15,159,
730
GMP_RNDD, "-2.5e09955c663d@65");
731
check2a(-1.49963910666191123860e+265,76,-2.30915090591874527520e-191,8,75,
732
GMP_RNDZ, "-1.dc3ec027da54e@220");
733
check2a(3.25471707846623300604e-160,81,-7.93846654265839958715e-274,58,54,
734
GMP_RNDN, "4.936a52bc17254@-133");
735
check2a(5.17945380930936917508e+112,119,1.11369077158813567738e+108,15,150,
736
GMP_RNDZ, "5.62661692498ec@93");
737
check2a(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,
738
GMP_RNDZ, "-a.204acdd25d788@-44");
739
check2a(-1.87427265794105342764e-57,175,1.76570844587489516446e+190,2,115,
740
GMP_RNDZ, "b.fffffffffffffffffffffffffffe@157");
741
check2a(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,
742
GMP_RNDU, "-b.eae2643497ff6286b@-108");
743
check2a(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,
744
GMP_RNDD, "-b.eae2643497ff6286b@-108");
745
check2a(-3.31624349995221499866e-22,107,-8.20150212714204839621e+156,79,99,
746
GMP_RNDD, "-2.63b22b55697e8000000000008@130");
747
x = -5943982715394951.0; for (i=0; i<446; i++) x *= 2.0;
748
check2a(x, 63, 1.77607317509426332389e+73, 64, 64, GMP_RNDN,
749
"-5.4781549356e1c@124");
750
check2a(4.49465557237618783128e+53,108,-2.45103927353799477871e+48,60,105,
751
GMP_RNDN, "4.b14f230f909dc803e@44");
752
check2a(2.26531902208967707071e+168,99,-2.67795218510613988524e+168,67,94,
753
GMP_RNDU, "-1.bfd7ff2647098@139");
754
check2a(-8.88471912490080158206e+253,79,-7.84488427404526918825e+124,95,53,
755
GMP_RNDD, "-c.1e533b8d835@210");
756
check2a(-2.18548638152863831959e-125,61,-1.22788940592412363564e+226,71,54,
757
GMP_RNDN, "-8.4b0f99ffa3b58@187");
758
check2a(-7.94156823309993162569e+77,74,-5.26820160805275124882e+80,99,101,
759
GMP_RNDD, "-1.1cc90f11d6af26f4@67");
760
check2a(-3.85170653452493859064e+189,62,2.18827389706660314090e+158,94,106,
761
GMP_RNDD, "-3.753ac0935b701ffffffffffffd@157");
762
check2a(1.07966151149311101950e+46,88,1.13198076934147457400e+46,67,53,
763
GMP_RNDN, "3.dfbc152dd4368@38");
764
check2a(3.36768223223409657622e+209,55,-9.61624007357265441884e+219,113,53,
765
GMP_RNDN, "-6.cf7217a451388@182");
766
check2a(-6.47376909368979326475e+159,111,5.11127211034490340501e+159,99,62,
767
GMP_RNDD, "-1.8cf3aadf537c@132");
768
check2a(-4.95229483271607845549e+220,110,-6.06992115033276044773e+213,109,55,
769
GMP_RNDN, "-2.3129f1f63b31b@183");
770
check2a(-6.47376909368979326475e+159,74,5.11127211034490340501e+159,111,75,
771
GMP_RNDU, "-1.8cf3aadf537c@132");
772
check2a(2.26531902208967707070e+168,99,-2.67795218510613988525e+168,67,94,
773
GMP_RNDU, "-1.bfd7ff2647098@139");
774
check2a(-2.28886326552077479586e-188,67,3.41419438647157839320e-177,60,110,
775
GMP_RNDU, "3.75740b4fe8f17f90258907@-147");
776
check2a(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,
777
GMP_RNDZ, "-a.204acdd25d788@-44");
778
check2a(2.90983392714730768886e+50,101,2.31299792168440591870e+50,74,105,
779
GMP_RNDZ, "1.655c53ff5719c8@42");
780
check2a(2.72046257722708717791e+243,97,-1.62158447436486437113e+243,83,96,
781
GMP_RNDN, "a.4cc63e002d2e8@201");
782
/* Checking double precision (53 bits) */
783
check53(-8.22183238641455905806e-19, 7.42227178769761587878e-19, GMP_RNDD,
784
-7.9956059871694317927e-20);
785
check53(5.82106394662028628236e+234, -5.21514064202368477230e+89, GMP_RNDD,
786
5.8210639466202855763e234);
787
check53(5.72931679569871602371e+122, -5.72886070363264321230e+122, GMP_RNDN,
788
4.5609206607281141508e118);
789
check53(-5.09937369394650450820e+238, 2.70203299854862982387e+250, GMP_RNDD,
790
2.7020329985435301323e250);
791
check53(-2.96695924472363684394e+27, 1.22842938251111500000e+16, GMP_RNDD,
792
-2.96695924471135255027e27);
793
check53(1.74693641655743793422e-227, -7.71776956366861843469e-229, GMP_RNDN,
794
1.669758720920751867e-227);
795
x = -7883040437021647.0; for (i=0; i<468; i++) x = x / 2.0;
796
check53(-1.03432206392780011159e-125, 1.30127034799251347548e-133, GMP_RNDN,
798
check53(1.05824655795525779205e+71, -1.06022698059744327881e+71, GMP_RNDZ,
799
-1.9804226421854867632e68);
800
check53(-5.84204911040921732219e+240, 7.26658169050749590763e+240, GMP_RNDD,
801
1.4245325800982785854e240);
802
check53(1.00944884131046636376e+221, 2.33809162651471520268e+215, GMP_RNDN,
803
1.0094511794020929787e221);
804
x = 7045852550057985.0; for (i=0; i<986; i++) x = x / 2.0;
805
check53(4.29232078932667367325e-278, x, GMP_RNDU,
806
4.2933981418314132787e-278);
807
check53(5.27584773801377058681e-80, 8.91207657803547196421e-91, GMP_RNDN,
808
5.2758477381028917269e-80);
809
check53(2.99280481918991653800e+272, 5.34637717585790933424e+271, GMP_RNDN,
810
3.5274425367757071711e272);
811
check53(4.67302514390488041733e-184, 2.18321376145645689945e-190, GMP_RNDN,
812
4.6730273271186420541e-184);
813
check53(5.57294120336300389254e+71, 2.60596167942024924040e+65, GMP_RNDZ,
814
5.5729438093246831053e71);
815
check53(6.6052588496951015469e24, 4938448004894539.0, GMP_RNDU,
816
6.6052588546335505068e24);
817
check53(1.23056185051606761523e-190, 1.64589756643433857138e-181, GMP_RNDU,
818
1.6458975676649006598e-181);
819
check53(2.93231171510175981584e-280, 3.26266919161341483877e-273, GMP_RNDU,
820
3.2626694848445867288e-273);
821
check53(5.76707395945001907217e-58, 4.74752971449827687074e-51, GMP_RNDD,
822
4.747530291205672325e-51);
823
check53(277363943109.0, 11.0, GMP_RNDN, 277363943120.0);
824
#if 0 /* disabled since it seems silly to use denorms *
825
/* test denormalized numbers too */
826
check53(8.06294740693074521573e-310, 6.95250701071929654575e-310, GMP_RNDU,
827
1.5015454417650041761e-309);
830
/* the following check double overflow */
831
check53(6.27557402141211962228e+307, 1.32141396570101687757e+308,
832
GMP_RNDZ, DBL_POS_INF);
833
check53(DBL_POS_INF, 6.95250701071929654575e-310, GMP_RNDU, DBL_POS_INF);
834
check53(DBL_NEG_INF, 6.95250701071929654575e-310, GMP_RNDU, DBL_NEG_INF);
835
check53(6.95250701071929654575e-310, DBL_POS_INF, GMP_RNDU, DBL_POS_INF);
836
check53(6.95250701071929654575e-310, DBL_NEG_INF, GMP_RNDU, DBL_NEG_INF);
837
check53nan (DBL_POS_INF, DBL_NEG_INF, GMP_RNDN);
839
check53(1.44791789689198883921e-140, -1.90982880222349071284e-121,
840
GMP_RNDN, -1.90982880222349071e-121);
843
/* tests for particular cases (Vincent Lefevre, 22 Aug 2001) */
844
check53(9007199254740992.0, 1.0, GMP_RNDN, 9007199254740992.0);
845
check53(9007199254740994.0, 1.0, GMP_RNDN, 9007199254740996.0);
846
check53(9007199254740992.0, -1.0, GMP_RNDN, 9007199254740991.0);
847
check53(9007199254740994.0, -1.0, GMP_RNDN, 9007199254740992.0);
848
check53(9007199254740996.0, -1.0, GMP_RNDN, 9007199254740996.0);
850
#ifdef MPFR_HAVE_FESETROUND
851
prec = (argc<2) ? 53 : atoi(argv[1]);
852
rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
853
/* Comparing to double precision using machine arithmetic */
857
if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
858
/* avoid denormalized numbers and overflows */
859
rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
860
check(x, y, rnd, prec, prec, prec, 0.0);
863
/* tests with random precisions */
866
px = 53 + (LONG_RAND() % 64);
867
py = 53 + (LONG_RAND() % 64);
868
pz = 53 + (LONG_RAND() % 64);
869
rnd_mode = LONG_RAND() % 4;
870
do { x = drand(); } while (isnan(x));
871
do { y = drand(); } while (isnan(y));
872
check2 (x, px, y, py, pz, rnd_mode);
874
/* Checking mpfr_add(x, x, y) with prec=53 */
878
if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
879
/* avoid denormalized numbers and overflows */
880
rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
884
/* Checking mpfr_add(x, y, x) with prec=53 */
888
if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
889
/* avoid denormalized numbers and overflows */
890
rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
894
/* Checking mpfr_add(x, x, x) with prec=53 */
896
do { x = drand(); } while ((ABS(x)<2.2e-307) || (ABS(x)>0.8e308));
897
/* avoid denormalized numbers and overflows */
898
rnd = (rnd_mode==-1) ? LONG_RAND()%4 : rnd_mode;
added
removed
src/gmp/mpfr/tests/tadd.c
