Annotation of OpenXM_contrib/gmp/mpfr/tests/tdiv.c, Revision 1.1.1.1
1.1 ohara 1: /* Test file for mpfr_div.
2:
3: Copyright 1999, 2001, 2002 Free Software Foundation, Inc.
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 <stdio.h>
23: #include <stdlib.h>
24: #include <time.h>
25: #include "gmp.h"
26: #include "gmp-impl.h"
27: #include "mpfr.h"
28: #include "mpfr-impl.h"
29: #include "mpfr-test.h"
30:
31: #define check53(n, d, rnd, res) check4(n, d, rnd, 53, res)
32:
33: void check4 _PROTO((double, double, mp_rnd_t, int, double));
34: void check24 _PROTO((float, float, mp_rnd_t, float));
35: void check_float _PROTO((void));
36: void check_convergence _PROTO((void));
37: void check_lowr _PROTO((void));
38: void check_inexact _PROTO((void));
39: void check_nan _PROTO((void));
40:
41: void
42: check4 (double N, double D, mp_rnd_t rnd_mode, int p, double Q)
43: {
44: mpfr_t q, n, d; double Q2;
45:
46: mpfr_init2(q, p); mpfr_init2(n, p); mpfr_init2(d, p);
47: mpfr_set_d(n, N, rnd_mode);
48: mpfr_set_d(d, D, rnd_mode);
49: mpfr_div(q, n, d, rnd_mode);
50: #ifdef MPFR_HAVE_FESETROUND
51: mpfr_set_machine_rnd_mode(rnd_mode);
52: #endif
53: if (Q==0.0) Q = N/D;
54: Q2 = mpfr_get_d1 (q);
55: if (p==53 && Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
56: printf("mpfr_div failed for n=%1.20e, d=%1.20e, rnd_mode=%s\n",
57: N, D, mpfr_print_rnd_mode(rnd_mode));
58: printf("expected quotient is %1.20e, got %1.20e (%d ulp)\n", Q, Q2,
59: ulp(Q2, Q));
60: exit(1);
61: }
62: mpfr_clear(q); mpfr_clear(n); mpfr_clear(d);
63: }
64:
65: void
66: check24 (float N, float D, mp_rnd_t rnd_mode, float Q)
67: {
68: mpfr_t q, n, d; float Q2;
69:
70: mpfr_init2(q, 24); mpfr_init2(n, 24); mpfr_init2(d, 24);
71: mpfr_set_d(n, N, rnd_mode);
72: mpfr_set_d(d, D, rnd_mode);
73: mpfr_div(q, n, d, rnd_mode);
74: Q2 = mpfr_get_d1 (q);
75: if (Q!=Q2) {
76: printf("mpfr_div failed for n=%1.10e, d=%1.10e, prec=24, rnd_mode=%s\n",
77: N, D, mpfr_print_rnd_mode(rnd_mode));
78: printf("expected quotient is %1.10e, got %1.10e\n", Q, Q2);
79: exit(1);
80: }
81: mpfr_clear(q); mpfr_clear(n); mpfr_clear(d);
82: }
83:
84: /* the following examples come from the paper "Number-theoretic Test
85: Generation for Directed Rounding" from Michael Parks, Table 2 */
86: void
87: check_float(void)
88: {
89: float b=8388608.0; /* 2^23 */
90:
91: check24(b*8388610.0, 8388609.0, GMP_RNDN, 8.388609e6);
92: check24(b*16777215.0, 16777213.0, GMP_RNDN, 8.388609e6);
93: check24(b*8388612.0, 8388611.0, GMP_RNDN, 8.388609e6);
94: check24(b*12582914.0, 12582911.0, GMP_RNDN, 8.38861e6);
95: check24(12582913.0, 12582910.0, GMP_RNDN, 1.000000238);
96: check24(b*16777215.0, 8388609.0, GMP_RNDN, 1.6777213e7);
97: check24(b*8388612.0, 8388609.0, GMP_RNDN, 8.388611e6);
98: check24(b*12582914.0, 8388610.0, GMP_RNDN, 1.2582911e7);
99: check24(b*12582913.0, 8388610.0, GMP_RNDN, 1.258291e7);
100:
101: check24(b*8388610.0, 8388609.0, GMP_RNDZ, 8.388608e6);
102: check24(b*16777215.0, 16777213.0, GMP_RNDZ, 8.388609e6);
103: check24(b*8388612.0, 8388611.0, GMP_RNDZ, 8.388608e6);
104: check24(b*12582914.0, 12582911.0, GMP_RNDZ, 8.38861e6);
105: check24(12582913.0, 12582910.0, GMP_RNDZ, 1.000000238);
106: check24(b*16777215.0, 8388609.0, GMP_RNDZ, 1.6777213e7);
107: check24(b*8388612.0, 8388609.0, GMP_RNDZ, 8.38861e6);
108: check24(b*12582914.0, 8388610.0, GMP_RNDZ, 1.2582911e7);
109: check24(b*12582913.0, 8388610.0, GMP_RNDZ, 1.258291e7);
110:
111: check24(b*8388610.0, 8388609.0, GMP_RNDU, 8.388609e6);
112: check24(b*16777215.0, 16777213.0, GMP_RNDU, 8.38861e6);
113: check24(b*8388612.0, 8388611.0, GMP_RNDU, 8.388609e6);
114: check24(b*12582914.0, 12582911.0, GMP_RNDU, 8.388611e6);
115: check24(12582913.0, 12582910.0, GMP_RNDU, 1.000000357);
116: check24(b*16777215.0, 8388609.0, GMP_RNDU, 1.6777214e7);
117: check24(b*8388612.0, 8388609.0, GMP_RNDU, 8.388611e6);
118: check24(b*12582914.0, 8388610.0, GMP_RNDU, 1.2582912e7);
119: check24(b*12582913.0, 8388610.0, GMP_RNDU, 1.2582911e7);
120:
121: check24(b*8388610.0, 8388609.0, GMP_RNDD, 8.388608e6);
122: check24(b*16777215.0, 16777213.0, GMP_RNDD, 8.388609e6);
123: check24(b*8388612.0, 8388611.0, GMP_RNDD, 8.388608e6);
124: check24(b*12582914.0, 12582911.0, GMP_RNDD, 8.38861e6);
125: check24(12582913.0, 12582910.0, GMP_RNDD, 1.000000238);
126: check24(b*16777215.0, 8388609.0, GMP_RNDD, 1.6777213e7);
127: check24(b*8388612.0, 8388609.0, GMP_RNDD, 8.38861e6);
128: check24(b*12582914.0, 8388610.0, GMP_RNDD, 1.2582911e7);
129: check24(b*12582913.0, 8388610.0, GMP_RNDD, 1.258291e7);
130: }
131:
132: void
133: check_convergence (void)
134: {
135: mpfr_t x, y; int i, j;
136:
137: mpfr_init2(x, 130);
138: mpfr_set_str_raw(x, "0.1011111101011010101000001010011111101000011100011101010011111011000011001010000000111100100111110011001010110100100001001000111001E6944");
139: mpfr_init2(y, 130);
140: mpfr_set_ui(y, 5, GMP_RNDN);
141: mpfr_div(x, x, y, GMP_RNDD); /* exact division */
142:
143: mpfr_set_prec(x, 64);
144: mpfr_set_prec(y, 64);
145: mpfr_set_str_raw(x, "0.10010010011011010100101001010111100000101110010010101E55");
146: mpfr_set_str_raw(y, "0.1E585");
147: mpfr_div(x, x, y, GMP_RNDN);
148: mpfr_set_str_raw(y, "0.10010010011011010100101001010111100000101110010010101E-529");
149: if (mpfr_cmp(x, y)) {
150: fprintf(stderr, "Error in mpfr_div for prec=64, rnd=GMP_RNDN\n");
151: printf("got "); mpfr_print_binary(x); putchar('\n');
152: printf("instead of "); mpfr_print_binary(y); putchar('\n');
153: exit(1);
154: }
155:
156: for (i=32; i<=64; i+=32) {
157: mpfr_set_prec(x, i);
158: mpfr_set_prec(y, i);
159: mpfr_set_d(x, 1.0, GMP_RNDN);
160: for (j=0;j<4;j++) {
161: mpfr_set_d(y, 1.0, GMP_RNDN);
162: mpfr_div(y, x, y, j);
163: if (mpfr_cmp_ui(y, 1)) {
164: fprintf(stderr, "mpfr_div failed for x=1.0, y=1.0, prec=%u rnd=%s\n",
165: i, mpfr_print_rnd_mode(j));
166: printf("got "); mpfr_print_binary(y); putchar('\n');
167: exit(1);
168: }
169: }
170: }
171:
172: mpfr_clear(x); mpfr_clear(y);
173: }
174:
175: void
176: check_lowr (void)
177: {
178: mpfr_t x, y, z, z2, z3, tmp;
179: int k, c;
180:
181:
182: mpfr_init2 (x, 1000);
183: mpfr_init2 (y, 100);
184: mpfr_init2 (tmp, 850);
185: mpfr_init2 (z, 10);
186: mpfr_init2 (z2, 10);
187: mpfr_init2 (z3, 50);
188:
189: for (k = 1; k < 10000; k++)
190: {
191: mpfr_random (z);
192: mpfr_random (tmp);
193: mpfr_mul (x, z, tmp, GMP_RNDN);
194: c = mpfr_div (z2, x, tmp, GMP_RNDN);
195:
196: if (c || mpfr_cmp(z2, z))
197: {
198: fprintf(stderr, "Error in mpfr_div rnd=GMP_RNDN\n");
199: printf("Dividing ");
200: printf("got "); mpfr_print_binary(z2); putchar('\n');
201: printf("instead of "); mpfr_print_binary(z); putchar('\n');
202: printf("inex flag = %d\n", c);
203: exit(1);
204: }
205: }
206:
207: mpfr_set_prec(z2, 9);
208: for (k = 1; k < 10000; k++)
209: {
210: mpfr_random(z);
211: mpfr_random(tmp);
212: mpfr_mul(x, z, tmp, GMP_RNDN);
213: c = mpfr_div(z2, x, tmp, GMP_RNDN);
214:
215: if ((mpfr_cmp(z2, z) == 0 && c) || c == -1)
216: {
217: fprintf(stderr, "Error in mpfr_div rnd=GMP_RNDN\n");
218: printf("Dividing ");
219: printf("got "); mpfr_print_binary(z2); putchar('\n');
220: printf("instead of "); mpfr_print_binary(z); putchar('\n');
221: printf("inex flag = %d\n", c);
222: exit(1);
223: }
224: else if (c == 2)
225: {
226: mpfr_add_one_ulp(z, GMP_RNDN);
227: if (mpfr_cmp(z2, z))
228: {
229: fprintf(stderr, "Error in mpfr_div [even rnd?] rnd=GMP_RNDN\n");
230: printf("Dividing ");
231: printf("got "); mpfr_print_binary(z2); putchar('\n');
232: printf("instead of "); mpfr_print_binary(z); putchar('\n');
233: printf("inex flag = %d\n", 1);
234: exit(1);
235: }
236: }
237: else if (c == -2)
238: {
239: mpfr_sub_one_ulp(z, GMP_RNDN);
240: if (mpfr_cmp(z2, z))
241: {
242: fprintf(stderr, "Error in mpfr_div [even rnd?] rnd=GMP_RNDN\n");
243: printf("Dividing ");
244: printf("got "); mpfr_print_binary(z2); putchar('\n');
245: printf("instead of "); mpfr_print_binary(z); putchar('\n');
246: printf("inex flag = %d\n", 1);
247: exit(1);
248: }
249: }
250: }
251:
252:
253: mpfr_set_prec(x, 1000);
254: mpfr_set_prec(y, 100);
255: mpfr_set_prec(tmp, 850);
256: mpfr_set_prec(z, 10);
257: mpfr_set_prec(z2, 10);
258:
259: /* almost exact divisions */
260: for (k = 1; k < 10000; k++)
261: {
262: mpfr_random(z);
263: mpfr_random(tmp);
264: mpfr_mul(x, z, tmp, GMP_RNDN);
265: mpfr_set(y, tmp, GMP_RNDD);
266: mpfr_add_one_ulp(x, GMP_RNDN);
267:
268: c = mpfr_div(z2, x, y, GMP_RNDD);
269: mpfr_div(z3, x, y, GMP_RNDD);
270: mpfr_set(z, z3, GMP_RNDD);
271:
272: if (c != -1 || mpfr_cmp(z2, z))
273: {
274: fprintf(stderr, "Error in mpfr_div rnd=GMP_RNDD\n");
275: printf("got "); mpfr_print_binary(z2); putchar('\n');
276: printf("instead of "); mpfr_print_binary(z); putchar('\n');
277: printf("inex flag = %d\n", c);
278: exit(1);
279: }
280:
281: mpfr_set(y, tmp, GMP_RNDU);
282: c = mpfr_div(z2, x, y, GMP_RNDU);
283: mpfr_div(z3, x, y, GMP_RNDU);
284: mpfr_set(z, z3, GMP_RNDU);
285: if (c != 1 || mpfr_cmp(z2, z))
286: {
287: fprintf(stderr, "Error in mpfr_div rnd=GMP_RNDU\n");
288: printf("got "); mpfr_print_binary(z2); putchar('\n');
289: printf("instead of "); mpfr_print_binary(z); putchar('\n');
290: printf("inex flag = %d\n", c);
291: exit(1);
292: }
293: }
294:
295: mpfr_clear (x);
296: mpfr_clear (y);
297: mpfr_clear (z);
298: mpfr_clear (z2);
299: mpfr_clear (z3);
300: mpfr_clear (tmp);
301: }
302:
303: #define MAX_PREC 100
304:
305: void
306: check_inexact (void)
307: {
308: mpfr_t x, y, z, u;
309: mp_prec_t px, py, pu;
310: int inexact, cmp;
311: mp_rnd_t rnd;
312:
313: mpfr_init (x);
314: mpfr_init (y);
315: mpfr_init (z);
316: mpfr_init (u);
317:
318: mpfr_set_prec (x, 33);
319: mpfr_set_str_raw (x, "0.101111100011011101010011101100001E0");
320: mpfr_set_prec (u, 2);
321: mpfr_set_str_raw (u, "0.1E0");
322: mpfr_set_prec (y, 28);
323: if ((inexact = mpfr_div (y, x, u, GMP_RNDN) >= 0))
324: {
325: fprintf (stderr, "Wrong inexact flag (1): expected -1, got %d\n",
326: inexact);
327: exit (1);
328: }
329:
330: mpfr_set_prec (x, 129);
331: mpfr_set_str_raw (x, "0.111110101111001100000101011100101100110011011101010001000110110101100101000010000001110110100001101010001010100010001111001101010E-2");
332: mpfr_set_prec (u, 15);
333: mpfr_set_str_raw (u, "0.101101000001100E-1");
334: mpfr_set_prec (y, 92);
335: if ((inexact = mpfr_div (y, x, u, GMP_RNDN) <= 0))
336: {
337: fprintf (stderr, "Wrong inexact flag (1): expected 1, got %d\n",
338: inexact);
339: mpfr_print_binary(y); putchar('\n');
340: exit (1);
341: }
342:
343: for (px=2; px<MAX_PREC; px++)
344: {
345: mpfr_set_prec (x, px);
346: mpfr_random (x);
347: for (pu=2; pu<=MAX_PREC; pu++)
348: {
349: mpfr_set_prec (u, pu);
350: do { mpfr_random (u); } while (mpfr_cmp_ui (u, 0) == 0);
351: for (py=2; py<=MAX_PREC; py++)
352: {
353: mpfr_set_prec (y, py);
354: mpfr_set_prec (z, py + pu);
355: for (rnd=0; rnd<4; rnd++)
356: {
357: inexact = mpfr_div (y, x, u, rnd);
358: if (mpfr_mul (z, y, u, rnd))
359: {
360: fprintf (stderr, "z <- y * u should be exact\n");
361: exit (1);
362: }
363: cmp = mpfr_cmp (z, x);
364: if (((inexact == 0) && (cmp != 0)) ||
365: ((inexact > 0) && (cmp <= 0)) ||
366: ((inexact < 0) && (cmp >= 0)))
367: {
368: fprintf (stderr, "Wrong inexact flag for rnd=%s\n",
369: mpfr_print_rnd_mode(rnd));
370: printf ("expected %d, got %d\n", cmp, inexact);
371: printf ("x="); mpfr_print_binary (x); putchar ('\n');
372: printf ("u="); mpfr_print_binary (u); putchar ('\n');
373: printf ("y="); mpfr_print_binary (y); putchar ('\n');
374: printf ("y*u="); mpfr_print_binary (z); putchar ('\n');
375: exit (1);
376: }
377: }
378: }
379: }
380: }
381:
382: mpfr_clear (x);
383: mpfr_clear (y);
384: mpfr_clear (z);
385: mpfr_clear (u);
386: }
387:
388: void
389: check_nan (void)
390: {
391: mpfr_t a, d, q;
392:
393: mpfr_init2 (a, 100L);
394: mpfr_init2 (d, 100L);
395: mpfr_init2 (q, 100L);
396:
397: /* 1/nan == nan */
398: mpfr_set_ui (a, 1L, GMP_RNDN);
399: MPFR_SET_NAN (d);
400: ASSERT_ALWAYS (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
401: ASSERT_ALWAYS (mpfr_nan_p (q));
402:
403: /* nan/1 == nan */
404: MPFR_SET_NAN (a);
405: mpfr_set_ui (d, 1L, GMP_RNDN);
406: ASSERT_ALWAYS (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
407: ASSERT_ALWAYS (mpfr_nan_p (q));
408:
409: /* +inf/1 == +inf */
410: MPFR_CLEAR_FLAGS (a);
411: MPFR_SET_INF (a);
412: MPFR_SET_POS (a);
413: mpfr_set_ui (d, 1L, GMP_RNDN);
414: ASSERT_ALWAYS (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
415: ASSERT_ALWAYS (mpfr_inf_p (q));
416: ASSERT_ALWAYS (mpfr_sgn (q) > 0);
417:
418: /* 1/+inf == 0 */
419: mpfr_set_ui (a, 1L, GMP_RNDN);
420: MPFR_CLEAR_FLAGS (d);
421: MPFR_SET_INF (d);
422: MPFR_SET_POS (d);
423: ASSERT_ALWAYS (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
424: ASSERT_ALWAYS (mpfr_number_p (q));
425: ASSERT_ALWAYS (mpfr_sgn (q) == 0);
426:
427: /* 0/0 == nan */
428: mpfr_set_ui (a, 0L, GMP_RNDN);
429: mpfr_set_ui (d, 0L, GMP_RNDN);
430: ASSERT_ALWAYS (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
431: ASSERT_ALWAYS (mpfr_nan_p (q));
432:
433: /* +inf/+inf == nan */
434: MPFR_CLEAR_FLAGS (a);
435: MPFR_SET_INF (a);
436: MPFR_SET_POS (a);
437: MPFR_CLEAR_FLAGS (d);
438: MPFR_SET_INF (d);
439: MPFR_SET_POS (d);
440: ASSERT_ALWAYS (mpfr_div (q, a, d, GMP_RNDZ) == 0); /* exact */
441: ASSERT_ALWAYS (mpfr_nan_p (q));
442:
443: mpfr_clear (a);
444: mpfr_clear (d);
445: mpfr_clear (q);
446: }
447:
448:
449: int
450: main (int argc, char *argv[])
451: {
452: mpfr_t x, y, z;
453:
454: #ifdef MPFR_HAVE_FESETROUND
455: int N, i;
456: double n, d, e;
457:
458: mpfr_test_init ();
459: #endif
460:
461: check_inexact();
462:
463: mpfr_init2 (x, 64);
464: mpfr_init2 (y, 64);
465: mpfr_init2 (z, 64);
466:
467: mpfr_set_str_raw(x, "1.00100100110110101001010010101111000001011100100101010000000000E54");
468: mpfr_set_str_raw(y, "1.00000000000000000000000000000000000000000000000000000000000000E584");
469: mpfr_div(z, x, y, GMP_RNDU);
470:
471: check_nan ();
472: check_lowr();
473: check_float(); /* checks single precision */
474: check_convergence();
475: check53(0.0, 1.0, GMP_RNDZ, 0.0);
476: check53(-7.4988969224688591e63, 4.8816866450288732e306, GMP_RNDD,
477: -1.5361282826510687291e-243);
478: check53(-1.33225773037748601769e+199, 3.63449540676937123913e+79, GMP_RNDZ,
479: -3.6655920045905428978e119);
480: check4(4.0, 4.503599627370496e15, GMP_RNDZ, 62, 0.0);
481: check4(1.0, 2.10263340267725788209e+187, GMP_RNDU, 65, 0.0);
482: check4(2.44394909079968374564e-150, 2.10263340267725788209e+187, GMP_RNDU,
483: 65, 0.0);
484: check53(9.89438396044940256501e-134, 5.93472984109987421717e-67, GMP_RNDU,
485: 1.6672003992376663654e-67);
486: check53(9.89438396044940256501e-134, -5.93472984109987421717e-67, GMP_RNDU,
487: -1.6672003992376663654e-67);
488: check53(-4.53063926135729747564e-308, 7.02293374921793516813e-84, GMP_RNDD,
489: -6.4512060388748850857e-225);
490: check53(6.25089225176473806123e-01, -2.35527154824420243364e-230, GMP_RNDD,
491: -2.6540006635008291192e229);
492: check53(6.52308934689126e15, -1.62063546601505417497e273, GMP_RNDN,
493: -4.0250194961676020848e-258);
494: check53(1.04636807108079349236e-189, 3.72295730823253012954e-292, GMP_RNDZ,
495: 2.810583051186143125e102);
496:
497: #ifdef MPFR_HAVE_FESETROUND
498: N = (argc>1) ? atoi(argv[1]) : 10000;
499: SEED_RAND (time(NULL));
500: for (i=0;i<N;i++)
501: {
502: do { n = drand(); d = drand(); e = ABS(n)/ABS(d); }
503: /* smallest normalized is 2^(-1022), largest is 2^(1023)*(2-2^(-52)) */
504: while (e>=MAXNORM || e<MINNORM);
505: check4 (n, d, LONG_RAND() % 4, 53, 0.0);
506: }
507: #endif
508:
509: mpfr_clear (x);
510: mpfr_clear (y);
511: mpfr_clear (z);
512:
513: return 0;
514: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tdiv.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