Annotation of OpenXM_contrib/gmp/mpfr/tests/tsub.c, Revision 1.1.1.1
1.1 ohara 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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tsub.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr / tests