Annotation of OpenXM_contrib/gmp/mpfr/tests/tpow3.c, Revision 1.1.1.1
1.1 ohara 1: /* Test file for mpfr_pow.
2:
3: Copyright 2001 Free Software Foundation.
4: Adapted from tarctan.c.
5:
6: This file is part of the MPFR Library.
7:
8: The MPFR Library is free software; you can redistribute it and/or modify
9: it under the terms of the GNU Lesser General Public License as published by
10: the Free Software Foundation; either version 2.1 of the License, or (at your
11: option) any later version.
12:
13: The MPFR Library is distributed in the hope that it will be useful, but
14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16: License for more details.
17:
18: You should have received a copy of the GNU Lesser General Public License
19: along with the MPFR Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: #include <stdio.h>
24: #include <limits.h>
25: #include <stdlib.h>
26: #include "gmp.h"
27: #include "gmp-impl.h"
28: #include "mpfr.h"
29: #include "mpfr-impl.h"
30: #include "mpfr-test.h"
31:
32:
33: int
34: main (int argc, char *argv[])
35: {
36: mpfr_t x, y, z, ax;
37: long int iy;
38: mpfr_init (x);
39: mpfr_init (ax);
40: mpfr_init2 (y,sizeof(unsigned long int)*CHAR_BIT);
41: mpfr_init (z);
42:
43: MPFR_SET_NAN(x);
44: mpfr_random(y);
45: mpfr_pow (z, x,y, GMP_RNDN);
46: if(!MPFR_IS_NAN(z))
47: {
48: printf ("evaluation of function in x=NAN does not return NAN");
49: exit (1);
50: }
51:
52: MPFR_SET_NAN(y);
53: mpfr_random(x);
54: mpfr_pow (z, x,y, GMP_RNDN);
55: if(!MPFR_IS_NAN(z))
56: {
57: printf ("evaluation of function in y=NAN does not return NAN");
58: exit (1);
59: }
60:
61: MPFR_CLEAR_FLAGS(z);
62: MPFR_CLEAR_FLAGS(y);
63: MPFR_CLEAR_FLAGS(x);
64:
65: MPFR_SET_ZERO(y);
66: mpfr_random(x);
67: mpfr_pow (z, x,y, GMP_RNDN);
68: if(mpfr_cmp_ui(z,1)!=0 && !(MPFR_IS_NAN(x)))
69: {
70: printf ("evaluation of function in y=0 does not return 1\n");
71: printf ("x =");
72: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
73: printf ("\n y =");
74: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
75: printf ("\n result =");
76: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
77: exit (1);
78: }
79:
80: MPFR_CLEAR_FLAGS(z);
81: MPFR_CLEAR_FLAGS(y);
82: MPFR_CLEAR_FLAGS(x);
83:
84: MPFR_SET_INF(y);
85: if (MPFR_SIGN(y) < 0)
86: MPFR_CHANGE_SIGN(y);
87: mpfr_random(x);
88: mpfr_set_prec (ax, MPFR_PREC(x));
89: mpfr_abs(ax,x,GMP_RNDN);
90: mpfr_pow (z, x,y, GMP_RNDN);
91: if( !MPFR_IS_INF(z) && (mpfr_cmp_ui(ax,1) > 0) )
92: {
93: printf ("evaluation of function in y=INF (|x|>1) does not return INF");
94: exit (1);
95: }
96: if( !MPFR_IS_ZERO(z) && (mpfr_cmp_ui(ax,1) < 0) )
97: {
98: printf ("\nevaluation of function in y=INF (|x|<1) does not return 0");
99: printf ("\nx =");
100: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
101: printf ("\n y =");
102: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
103: printf ("\n result =");
104: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
105: putchar('\n');
106: exit (1);
107: }
108:
109:
110: MPFR_CLEAR_FLAGS(z);
111: MPFR_CLEAR_FLAGS(y);
112: MPFR_CLEAR_FLAGS(x);
113:
114: MPFR_SET_INF(y);
115: if (MPFR_SIGN(y) > 0)
116: MPFR_CHANGE_SIGN(y);
117: mpfr_random(x);
118: mpfr_set_prec (ax, MPFR_PREC(x));
119: mpfr_abs(ax,x,GMP_RNDN);
120: mpfr_pow (z, x,y, GMP_RNDN);
121: mpfr_pow (z, x,y, GMP_RNDN);
122: if( !MPFR_IS_INF(z) && (mpfr_cmp_ui(ax,1) < 0) )
123: {
124: printf ("evaluation of function in y=INF (for |x| <0) does not return INF");
125: exit (1);
126: }
127: if( !MPFR_IS_ZERO(z) && (mpfr_cmp_ui(ax,1) > 0) )
128: {
129: printf ("evaluation of function in y=INF (for |x| >0) does not return 0");
130: exit (1);
131: }
132:
133: MPFR_CLEAR_FLAGS(z);
134: MPFR_CLEAR_FLAGS(y);
135: MPFR_CLEAR_FLAGS(x);
136:
137: MPFR_SET_INF(x);
138: if (MPFR_SIGN(x) < 0)
139: MPFR_CHANGE_SIGN(x);
140: mpfr_random(y);
141: mpfr_pow (z, x,y, GMP_RNDN);
142: if(!MPFR_IS_INF(z) && (MPFR_SIGN(y) > 0))
143: {
144: printf ("evaluation of function in INF does not return INF");
145: printf ("\nx =");
146: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
147: printf ("\n y =");
148: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
149: printf ("\n result =");
150: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
151: putchar('\n');
152: exit (1);
153: }
154: if(!MPFR_IS_ZERO(z) && (MPFR_SIGN(y) < 0))
155: {
156: printf ("evaluation of function in INF does not return INF");
157: printf ("\nx =");
158: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
159: printf ("\n y =");
160: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
161: printf ("\n result =");
162: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
163: putchar('\n');
164: exit (1);
165: }
166:
167:
168: MPFR_CLEAR_FLAGS(z);
169: MPFR_CLEAR_FLAGS(y);
170: MPFR_CLEAR_FLAGS(x);
171:
172: MPFR_SET_INF(x);
173: if (MPFR_SIGN(x) > 0)
174: MPFR_CHANGE_SIGN(x);
175: mpfr_random(y);
176: if (random() % 2)
177: mpfr_neg (y, y, GMP_RNDN);
178: mpfr_pow (z, x,y, GMP_RNDN);
179: if(!MPFR_IS_INF(z) && (MPFR_SIGN(y) > 0) && (mpfr_isinteger(y)))
180: {
181: printf ("evaluation of function in x=-INF does not return INF");
182: printf ("\nx =");
183: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
184: printf ("\n y =");
185: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
186: printf ("\n result =");
187: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
188: putchar('\n');
189: if(mpfr_isinteger(y))
190: printf("y is an integer\n");
191: else
192: printf("y is not an integer\n");
193:
194: exit (1);
195: }
196: if(!MPFR_IS_ZERO(z) && (MPFR_SIGN(y) < 0) && (mpfr_isinteger(y)))
197: {
198: printf ("evaluation of function in x=-INF does not return 0");
199: printf ("\nx =");
200: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
201: printf ("\n y =");
202: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
203: printf ("\n result =");
204: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
205: putchar('\n');
206:
207: if(mpfr_isinteger(y))
208: printf("y is an integer\n");
209: else
210: printf("y is not an integer\n");
211:
212: exit (1);
213: }
214: MPFR_CLEAR_FLAGS(z);
215: MPFR_CLEAR_FLAGS(y);
216: MPFR_CLEAR_FLAGS(x);
217:
218: MPFR_SET_INF(x);
219: if (MPFR_SIGN(x) > 0)
220: MPFR_CHANGE_SIGN(x);
221:
222: iy=random();
223: mpfr_random(y);
224: if (random() % 2)
225: iy=-iy;
226: mpfr_set_d(y,iy,GMP_RNDN);
227: mpfr_pow (z, x,y, GMP_RNDN);
228: if(!MPFR_IS_INF(z) && (MPFR_SIGN(y) > 0) && (mpfr_isinteger(y)))
229: {
230: printf ("evaluation of function in x=-INF does not return INF");
231: printf ("\nx =");
232: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
233: printf ("\n y =");
234: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
235: printf ("\n result =");
236: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
237: putchar('\n');
238: if(mpfr_isinteger(y))
239: printf("y is an integer\n");
240: else
241: printf("y is not an integer\n");
242:
243: exit (1);
244: }
245: if(!MPFR_IS_ZERO(z) && (MPFR_SIGN(y) < 0) && (mpfr_isinteger(y)))
246: {
247: printf ("evaluation of function in x=-INF does not return 0");
248: printf ("\nx =");
249: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
250: printf ("\n y =");
251: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
252: printf ("\n result =");
253: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
254: putchar('\n');
255:
256: if(mpfr_isinteger(y))
257: printf("y is an integer\n");
258: else
259: printf("y is not an integer\n");
260:
261: exit (1);
262: }
263: MPFR_CLEAR_FLAGS(z);
264: MPFR_CLEAR_FLAGS(y);
265: MPFR_CLEAR_FLAGS(x);
266:
267: mpfr_set_ui(x,1,GMP_RNDN);
268: MPFR_SET_INF(y);
269: mpfr_pow (z, x,y, GMP_RNDN);
270: if(!MPFR_IS_NAN(z))
271: {
272: printf ("evaluation of function in x=1, y=INF does not return NAN");
273: printf ("\nx =");
274: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
275: printf ("\n y =");
276: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
277: printf ("\n result =");
278: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
279: putchar('\n');
280:
281: exit (1);
282: }
283:
284: MPFR_CLEAR_FLAGS(z);
285: MPFR_CLEAR_FLAGS(y);
286: MPFR_CLEAR_FLAGS(x);
287:
288: MPFR_SET_ZERO(x);
289: mpfr_random(y);
290: if (random() % 2)
291: mpfr_neg (y, y, GMP_RNDN);
292:
293: mpfr_pow (z, x,y, GMP_RNDN);
294:
295: if(!MPFR_IS_ZERO(z) && (MPFR_SIGN(y) < 0) && !(mpfr_isinteger(y)))
296: {
297: printf ("evaluation of function in y<0 does not return 0");
298: printf ("\nx =");
299: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
300: printf ("\n y =");
301: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
302: printf ("\n result =");
303: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
304: putchar('\n');
305:
306: exit (1);
307: }
308: if(!MPFR_IS_INF(z) && (MPFR_SIGN(y) < 0) && (mpfr_isinteger(y)))
309: {
310: printf ("evaluation of function in y<0 (y integer) does not return INF");
311: printf ("\nx =");
312: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
313: printf ("\n y =");
314: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
315: printf ("\n result =");
316: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
317: putchar('\n');
318: exit (1);
319: }
320: if(!MPFR_IS_ZERO(z) && (MPFR_SIGN(y) > 0) && (mpfr_isinteger(y)))
321: {
322: printf ("evaluation of function in y<0 (y integer) does not return 0");
323: printf ("\nx =");
324: mpfr_out_str (stdout, 10, MPFR_PREC(x), x, GMP_RNDN);
325: printf ("\n y =");
326: mpfr_out_str (stdout, 10, MPFR_PREC(y), y, GMP_RNDN);
327: printf ("\n result =");
328: mpfr_out_str (stdout, 10, MPFR_PREC(z), z, GMP_RNDN);
329: putchar('\n');
330: exit (1);
331: }
332:
333:
334: {
335: mp_prec_t prec, yprec;
336: mpfr_t t, s;
337: mp_rnd_t rnd;
338: int inexact, compare, compare2;
339: unsigned int n, err;
340:
341: int p0=2;
342: int p1=100;
343: int N=100;
344:
345: mpfr_init (s);
346: mpfr_init (t);
347:
348: /* generic test */
349: for (prec = p0; prec <= p1; prec++)
350: {
351: mpfr_set_prec (x, prec);
352: mpfr_set_prec (s, sizeof(unsigned long int)*CHAR_BIT);
353: mpfr_set_prec (z, prec);
354: mpfr_set_prec (t, prec);
355: yprec = prec + 10;
356:
357: for (n=0; n<N; n++)
358: {
359:
360: mpfr_random (x);
361:
362: mpfr_random (s);
363: if (random() % 2)
364: mpfr_neg (s, s, GMP_RNDN);
365: rnd = random () % 4;
366: mpfr_set_prec (y, yprec);
367: compare = mpfr_pow (y, x, s, rnd);
368: err = (rnd == GMP_RNDN) ? yprec + 1 : yprec;
369: if (mpfr_can_round (y, err, rnd, rnd, prec))
370: {
371: mpfr_set (t, y, rnd);
372: inexact = mpfr_pow (z,x, s, rnd);
373: if (mpfr_cmp (t, z))
374: {
375: printf ("results differ for x=");
376: mpfr_out_str (stdout, 2, prec, x, GMP_RNDN);
377: printf (" values of the exponential=");
378: mpfr_out_str (stdout, 2, prec, s, GMP_RNDN);
379: printf (" prec=%u rnd_mode=%s\n", (unsigned) prec,
380: mpfr_print_rnd_mode (rnd));
381: printf ("got ");
382: mpfr_out_str (stdout, 2, prec, z, GMP_RNDN);
383: putchar ('\n');
384: printf ("expected ");
385: mpfr_out_str (stdout, 2, prec, t, GMP_RNDN);
386: putchar ('\n');
387: printf ("approx ");
388: mpfr_print_binary (y);
389: putchar ('\n');
390: exit (1);
391: }
392: compare2 = mpfr_cmp (t, y);
393: /* if rounding to nearest, cannot know the sign of t - f(x)
394: because of composed rounding: y = o(f(x)) and t = o(y) */
395: if ((rnd != GMP_RNDN) && (compare * compare2 >= 0))
396: compare = compare + compare2;
397: else
398: compare = inexact; /* cannot determine sign(t-f(x)) */
399: if (((inexact == 0) && (compare != 0)) ||
400: ((inexact > 0) && (compare <= 0)) ||
401: ((inexact < 0) && (compare >= 0)))
402: {
403: fprintf (stderr, "Wrong inexact flag for rnd=%s: expected %d, got %d\n",
404: mpfr_print_rnd_mode (rnd), compare, inexact);
405: printf ("x="); mpfr_print_binary (x); putchar ('\n');
406: printf ("y="); mpfr_print_binary (y); putchar ('\n');
407: printf ("t="); mpfr_print_binary (t); putchar ('\n');
408: exit (1);
409: }
410: }
411: }
412: }
413:
414: mpfr_clear (s);
415: mpfr_clear (t);
416:
417: }
418: mpfr_clear (x);
419: mpfr_clear (y);
420: mpfr_clear (z);
421: mpfr_clear (ax);
422:
423: return 0;
424: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>