Annotation of OpenXM_contrib/gmp/mpfr/tests/tsqrt.c, Revision 1.1.1.1
1.1 ohara 1: /* Test file for mpfr_sqrt.
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 <math.h>
23: #include <stdio.h>
24: #include <stdlib.h>
25: #include <time.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: #define check(a,r) check3(a,r,-1.0)
33:
34: int maxulp=0;
35:
36: void check3 _PROTO((double, mp_rnd_t, double));
37: void check4 _PROTO((double, mp_rnd_t, char *));
38: void check24 _PROTO((float, mp_rnd_t, float));
39: void check_float _PROTO((void));
40: void special _PROTO((void));
41: void check_inexact _PROTO((mp_prec_t));
42: void check_nan _PROTO((void));
43:
44: void
45: check3 (double a, mp_rnd_t rnd_mode, double Q)
46: {
47: mpfr_t q; double Q2; int ck,u;
48:
49: ck = (Q!=-1.0); /* if ck=1, then Q is certified correct */
50: mpfr_init2(q, 53);
51: mpfr_set_d(q, a, rnd_mode);
52: #ifdef MPFR_HAVE_FESETROUND
53: mpfr_set_machine_rnd_mode(rnd_mode);
54: #endif
55: mpfr_sqrt(q, q, rnd_mode);
56: if (ck==0) Q = sqrt(a);
57: else {
58: if (Q != sqrt(a) && (!isnan(Q) || !isnan(sqrt(a)))) {
59: fprintf(stderr, "you've found a bug in your machine's sqrt for x=%1.20e\n", a);
60: mpfr_clear(q);
61: exit(1);
62:
63: }
64: }
65: Q2 = mpfr_get_d1 (q);
66: if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
67: u = ulp(Q2,Q);
68: if (ck) {
69: printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%s\n",
70: a, mpfr_print_rnd_mode(rnd_mode));
71: printf("expected sqrt is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,u);
72: mpfr_clear(q);
73: exit(1);
74: }
75: else if (u>maxulp || u<-maxulp) {
76: maxulp = (u>maxulp) ? u : -u;
77: printf("libm.a differs from mpfr_sqrt for a=%1.20e, rnd_mode=%s\n",
78: a, mpfr_print_rnd_mode(rnd_mode));
79: printf("libm.a gives %1.20e, mpfr_sqrt gives %1.20e (%d ulp)\n",Q,Q2,u);
80: }
81: }
82: mpfr_clear(q);
83: }
84:
85: void
86: check4 (double a, mp_rnd_t rnd_mode, char *Q)
87: {
88: mpfr_t q, res;
89:
90: mpfr_init2(q, 53); mpfr_init2(res, 53);
91: mpfr_set_d(q, a, rnd_mode);
92: mpfr_sqrt(q, q, rnd_mode);
93: mpfr_set_str(res, Q, 16, GMP_RNDN);
94: if (mpfr_cmp(q, res)) {
95: printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%s\n",
96: a, mpfr_print_rnd_mode(rnd_mode));
97: printf("expected "); mpfr_print_binary(res); putchar('\n');
98: printf("got "); mpfr_print_binary(q); putchar('\n');
99: mpfr_clear(q); mpfr_clear(res);
100: exit(1);
101: }
102: mpfr_clear(res);
103: mpfr_clear(q);
104: }
105:
106: void
107: check24 (float a, mp_rnd_t rnd_mode, float Q)
108: {
109: mpfr_t q; float Q2;
110:
111: mpfr_init2(q, 24);
112: mpfr_set_d(q, a, rnd_mode);
113: mpfr_sqrt(q, q, rnd_mode);
114: Q2 = mpfr_get_d1 (q);
115: if (Q!=Q2) {
116: printf("mpfr_sqrt failed for a=%1.10e, prec=24, rnd_mode=%s\n",
117: a, mpfr_print_rnd_mode(rnd_mode));
118: printf("expected sqrt is %1.10e, got %1.10e\n",Q,Q2);
119: exit(1);
120: }
121: mpfr_clear(q);
122: }
123:
124: /* the following examples come from the paper "Number-theoretic Test
125: Generation for Directed Rounding" from Michael Parks, Table 3 */
126: void
127: check_float (void)
128: {
129: float b = 8388608.0; /* 2^23 */
130:
131: check24(b*8388610.0, GMP_RNDN, 8.388609e6);
132: check24(b*2.0*16777214.0, GMP_RNDN, 1.6777215e7);
133: check24(b*8388612.0, GMP_RNDN, 8.388610e6);
134: check24(b*2.0*16777212.0, GMP_RNDN, 1.6777214e7);
135: check24(b*11946704.0, GMP_RNDN, 1.0010805e7);
136: check24(b*14321479.0, GMP_RNDN, 1.0960715e7);
137: check24(b*2.0*13689673.0, GMP_RNDN, 1.5155019e7);
138: check24(b*8388614.0, GMP_RNDN, 8.388611e6);
139: check24(b*2.0*16777210.0, GMP_RNDN, 1.6777213e7);
140: check24(b*10873622.0, GMP_RNDN, 9.550631e6);
141:
142: check24(b*8388610.0, GMP_RNDZ, 8.388608e6);
143: check24(b*2.0*16777214.0, GMP_RNDZ, 1.6777214e7);
144: check24(b*8388612.0, GMP_RNDZ, 8.388609e6);
145: check24(b*2.0*16777212.0, GMP_RNDZ, 1.6777213e7);
146: check24(b*11946704.0, GMP_RNDZ, 1.0010805e7);
147: check24(b*14321479.0, GMP_RNDZ, 1.0960715e7);
148: check24(b*2.0*13689673.0, GMP_RNDZ, 1.5155019e7);
149: check24(b*8388614.0, GMP_RNDZ, 8.38861e6);
150: check24(b*2.0*16777210.0, GMP_RNDZ, 1.6777212e7);
151: check24(b*10873622.0, GMP_RNDZ, 9.550631e6);
152:
153: check24(b*8388610.0, GMP_RNDU, 8.388609e6);
154: check24(b*2.0*16777214.0, GMP_RNDU, 1.6777215e7);
155: check24(b*8388612.0, GMP_RNDU, 8.388610e6);
156: check24(b*2.0*16777212.0, GMP_RNDU, 1.6777214e7);
157: check24(b*11946704.0, GMP_RNDU, 1.0010806e7);
158: check24(b*14321479.0, GMP_RNDU, 1.0960716e7);
159: check24(b*2.0*13689673.0, GMP_RNDU, 1.515502e7);
160: check24(b*8388614.0, GMP_RNDU, 8.388611e6);
161: check24(b*2.0*16777210.0, GMP_RNDU, 1.6777213e7);
162: check24(b*10873622.0, GMP_RNDU, 9.550632e6);
163:
164: check24(b*8388610.0, GMP_RNDD, 8.388608e6);
165: check24(b*2.0*16777214.0, GMP_RNDD, 1.6777214e7);
166: check24(b*8388612.0, GMP_RNDD, 8.388609e6);
167: check24(b*2.0*16777212.0, GMP_RNDD, 1.6777213e7);
168: check24(b*11946704.0, GMP_RNDD, 1.0010805e7);
169: check24(b*14321479.0, GMP_RNDD, 1.0960715e7);
170: check24(b*2.0*13689673.0, GMP_RNDD, 1.5155019e7);
171: check24(b*8388614.0, GMP_RNDD, 8.38861e6);
172: check24(b*2.0*16777210.0, GMP_RNDD, 1.6777212e7);
173: check24(b*10873622.0, GMP_RNDD, 9.550631e6);
174: }
175:
176: void
177: special (void)
178: {
179: mpfr_t x, z;
180: int inexact;
181: mp_prec_t p;
182:
183: mpfr_init (x);
184: mpfr_init (z);
185:
186: mpfr_set_prec (x, 27);
187: mpfr_set_str_raw (x, "0.110100111010101000010001011");
188: if ((inexact = mpfr_sqrt (x, x, GMP_RNDZ)) >= 0)
189: {
190: fprintf (stderr, "Wrong inexact flag: expected -1, got %d\n", inexact);
191: exit (1);
192: }
193:
194: mpfr_set_prec (x, 2);
195: for (p=2; p<1000; p++)
196: {
197: mpfr_set_prec (z, p);
198: mpfr_set_ui (z, 1, GMP_RNDN);
199: mpfr_add_one_ulp (z, GMP_RNDN);
200: mpfr_sqrt (x, z, GMP_RNDU);
201: if (mpfr_get_d1 (x) != 1.5)
202: {
203: fprintf (stderr, "Error: sqrt(1+ulp(1), up) should give 1.5 (prec=%u)\n", (unsigned) p);
204: printf ("got "); mpfr_print_binary (x); putchar ('\n');
205: exit (1);
206: }
207: }
208:
209: /* check inexact flag */
210: mpfr_set_prec (x, 5);
211: mpfr_set_str_raw (x, "1.1001E-2");
212: if ((inexact = mpfr_sqrt (x, x, GMP_RNDN)))
213: {
214: fprintf (stderr, "Wrong inexact flag: expected 0, got %d\n", inexact);
215: exit (1);
216: }
217:
218: mpfr_set_prec (x, 2);
219: mpfr_set_prec (z, 2);
220:
221: /* checks the sign is correctly set */
222: mpfr_set_d (x, 1.0, GMP_RNDN);
223: mpfr_set_d (z, -1.0, GMP_RNDN);
224: mpfr_sqrt (z, x, GMP_RNDN);
225: if (mpfr_cmp_ui (z, 0) < 0) {
226: fprintf (stderr, "Error: square root of %e gives %e\n",
227: mpfr_get_d1 (x), mpfr_get_d1 (z));
228: exit (1);
229: }
230:
231: mpfr_set_prec (x, 192);
232: mpfr_set_prec (z, 160);
233: mpfr_set_str_raw (z, "0.1011010100000100100100100110011001011100100100000011000111011001011101101101110000110100001000100001100001011000E1");
234: mpfr_set_prec (x, 160);
235: mpfr_sqrt(x, z, GMP_RNDN);
236: mpfr_sqrt(z, x, GMP_RNDN);
237:
238: mpfr_clear (x);
239: mpfr_clear (z);
240: }
241:
242: void
243: check_inexact (mp_prec_t p)
244: {
245: mpfr_t x, y, z;
246: mp_rnd_t rnd;
247: int inexact, sign;
248:
249: mpfr_init2 (x, p);
250: mpfr_init2 (y, p);
251: mpfr_init2 (z, 2*p);
252: mpfr_random (x);
253: rnd = LONG_RAND() % 4;
254: inexact = mpfr_sqrt (y, x, rnd);
255: if (mpfr_mul (z, y, y, rnd)) /* exact since prec(z) = 2*prec(y) */
256: {
257: fprintf (stderr, "Error: multiplication should be exact\n");
258: exit (1);
259: }
260: mpfr_sub (z, z, x, rnd); /* exact also */
261: sign = mpfr_cmp_ui (z, 0);
262: if (((inexact == 0) && (sign)) ||
263: ((inexact > 0) && (sign <= 0)) ||
264: ((inexact < 0) && (sign >= 0)))
265: {
266: fprintf (stderr, "Error: wrong inexact flag, expected %d, got %d\n",
267: sign, inexact);
268: printf ("x=");
269: mpfr_print_binary (x);
270: printf (" rnd=%s\n", mpfr_print_rnd_mode (rnd));
271: printf ("y="); mpfr_print_binary (y); putchar ('\n');
272: exit (1);
273: }
274: mpfr_clear (x);
275: mpfr_clear (y);
276: mpfr_clear (z);
277: }
278:
279: void
280: check_nan (void)
281: {
282: mpfr_t x, got;
283:
284: mpfr_init2 (x, 100L);
285: mpfr_init2 (got, 100L);
286:
287: /* sqrt(NaN) == NaN */
288: MPFR_CLEAR_FLAGS (x);
289: MPFR_SET_NAN (x);
290: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
291: ASSERT_ALWAYS (mpfr_nan_p (got));
292:
293: /* sqrt(-1) == NaN */
294: mpfr_set_si (x, -1L, GMP_RNDZ);
295: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
296: ASSERT_ALWAYS (mpfr_nan_p (got));
297:
298: /* sqrt(+inf) == +inf */
299: MPFR_CLEAR_FLAGS (x);
300: MPFR_SET_INF (x);
301: MPFR_SET_POS (x);
302: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
303: ASSERT_ALWAYS (mpfr_inf_p (got));
304:
305: /* sqrt(-inf) == NaN */
306: MPFR_CLEAR_FLAGS (x);
307: MPFR_SET_INF (x);
308: MPFR_SET_NEG (x);
309: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
310: ASSERT_ALWAYS (mpfr_nan_p (got));
311:
312: /* sqrt(-0) == 0 */
313: mpfr_set_si (x, 0L, GMP_RNDZ);
314: MPFR_SET_NEG (x);
315: ASSERT_ALWAYS (mpfr_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
316: ASSERT_ALWAYS (mpfr_number_p (got));
317: ASSERT_ALWAYS (mpfr_cmp_ui (got, 0L) == 0);
318:
319: mpfr_clear (x);
320: mpfr_clear (got);
321: }
322:
323: double five = 5.0;
324:
325: int
326: main (void)
327: {
328: double a;
329: mp_prec_t p;
330: int k;
331: #ifdef MPFR_HAVE_FESETROUND
332: int i;
333:
334: mpfr_test_init ();
335:
336: /* On Debian potato glibc 2.1.3-18, sqrt() doesn't seem to respect
337: fesetround. */
338: {
339: double a, b;
340: mpfr_set_machine_rnd_mode (GMP_RNDU);
341: a = sqrt (five);
342: mpfr_set_machine_rnd_mode (GMP_RNDD);
343: b = sqrt (five);
344: if (a == b)
345: {
346: printf ("Tests suppressed, mpfr_set_machine_rnd_mode doesn't affect sqrt()\n");
347: goto nogood;
348: }
349: }
350:
351: SEED_RAND (time(NULL));
352: for (i=0;i<100000;i++)
353: {
354: a = drand();
355: if (a < 0.0) a = -a; /* ensures a is positive */
356: check (a, LONG_RAND() % 4);
357: }
358: nogood:
359: #endif
360:
361: check_nan ();
362:
363: for (p=2; p<200; p++)
364: for (k=0; k<200; k++)
365: check_inexact (p);
366: special ();
367: check_float();
368: #ifdef HAVE_INFS
369: check3 (DBL_NAN, GMP_RNDN, DBL_NAN);
370: check3 (-1.0, GMP_RNDN, DBL_NAN);
371: check3 (DBL_POS_INF, GMP_RNDN, DBL_POS_INF);
372: check3 (DBL_NEG_INF, GMP_RNDN, DBL_NAN);
373: #endif
374: check3(-0.0, GMP_RNDN, 0.0);
375: check4(6.37983013646045901440e+32, GMP_RNDN, "5.9bc5036d09e0c@13");
376: check4(1.0, GMP_RNDN, "1");
377: check4(1.0, GMP_RNDZ, "1");
378: check4(3.725290298461914062500000e-9, GMP_RNDN, "4@-4");
379: check4(3.725290298461914062500000e-9, GMP_RNDZ, "4@-4");
380: a=1190456976439861.0;
381: check4(a, GMP_RNDZ, "2.0e7957873529a@6");
382: check4(1024.0*a, GMP_RNDZ, "4.1cf2af0e6a534@7");
383: /* the following examples are bugs in Cygnus compiler/system, found by
384: Fabrice Rouillier while porting mpfr to Windows */
385: check4(9.89438396044940256501e-134, GMP_RNDU, "8.7af7bf0ebbee@-56");
386: check4(7.86528588050363751914e+31, GMP_RNDZ, "1.f81fc40f32062@13");
387: check4(0.99999999999999988897, GMP_RNDN, "f.ffffffffffff8@-1");
388: check4(1.00000000000000022204, GMP_RNDN, "1");
389: /* the following examples come from the paper "Number-theoretic Test
390: Generation for Directed Rounding" from Michael Parks, Table 4 */
391: a = 4503599627370496.0; /* 2^52 */
392:
393: check4(a*2.0*8732221479794286.0, GMP_RNDN, "1.f81fc40f32063@13");
394: check4(a*8550954388695124.0, GMP_RNDN, "1.60c012a92fc65@13");
395: check4(a*7842344481681754.0, GMP_RNDN, "1.51d17526c7161@13");
396: check4(a*5935035262218600.0, GMP_RNDN, "1.25e19302f7e51@13");
397: check4(a*5039650445085418.0, GMP_RNDN, "1.0ecea7dd2ec3d@13");
398: check4(a*5039721545366078.0, GMP_RNDN, "1.0ecf250e8e921@13");
399: check4(a*8005963117781324.0, GMP_RNDN, "1.5552f3eedcf33@13");
400: check4(a*6703494707970582.0, GMP_RNDN, "1.3853ee10c9c99@13");
401: check4(a*8010323124937260.0, GMP_RNDN, "1.556abe212b56f@13");
402: check4(a*2.0*8010776873384260.0, GMP_RNDN, "1.e2d9a51977e6e@13");
403:
404: check4(a*2.0*8732221479794286.0, GMP_RNDZ, "1.f81fc40f32062@13");
405: check4(a*8550954388695124.0, GMP_RNDZ, "1.60c012a92fc64@13");
406: check4(a*7842344481681754.0, GMP_RNDZ, "1.51d17526c716@13");
407: check4(a*5935035262218600.0, GMP_RNDZ, "1.25e19302f7e5@13");
408: check4(a*5039650445085418.0, GMP_RNDZ, "1.0ecea7dd2ec3c@13");
409: check4(a*5039721545366078.0, GMP_RNDZ, "1.0ecf250e8e92@13");
410: check4(a*8005963117781324.0, GMP_RNDZ, "1.5552f3eedcf32@13");
411: check4(a*6703494707970582.0, GMP_RNDZ, "1.3853ee10c9c98@13");
412: check4(a*8010323124937260.0, GMP_RNDZ, "1.556abe212b56e@13");
413: check4(a*2.0*8010776873384260.0, GMP_RNDZ, "1.e2d9a51977e6d@13");
414:
415: check4(a*2.0*8732221479794286.0, GMP_RNDU, "1.f81fc40f32063@13");
416: check4(a*8550954388695124.0, GMP_RNDU, "1.60c012a92fc65@13");
417: check4(a*7842344481681754.0, GMP_RNDU, "1.51d17526c7161@13");
418: check4(a*5935035262218600.0, GMP_RNDU, "1.25e19302f7e51@13");
419: check4(a*5039650445085418.0, GMP_RNDU, "1.0ecea7dd2ec3d@13");
420: check4(a*5039721545366078.0, GMP_RNDU, "1.0ecf250e8e921@13");
421: check4(a*8005963117781324.0, GMP_RNDU, "1.5552f3eedcf33@13");
422: check4(a*6703494707970582.0, GMP_RNDU, "1.3853ee10c9c99@13");
423: check4(a*8010323124937260.0, GMP_RNDU, "1.556abe212b56f@13");
424: check4(a*2.0*8010776873384260.0, GMP_RNDU, "1.e2d9a51977e6e@13");
425:
426: check4(a*2.0*8732221479794286.0, GMP_RNDD, "1.f81fc40f32062@13");
427: check4(a*8550954388695124.0, GMP_RNDD, "1.60c012a92fc64@13");
428: check4(a*7842344481681754.0, GMP_RNDD, "1.51d17526c716@13");
429: check4(a*5935035262218600.0, GMP_RNDD, "1.25e19302f7e5@13");
430: check4(a*5039650445085418.0, GMP_RNDD, "1.0ecea7dd2ec3c@13");
431: check4(a*5039721545366078.0, GMP_RNDD, "1.0ecf250e8e92@13");
432: check4(a*8005963117781324.0, GMP_RNDD, "1.5552f3eedcf32@13");
433: check4(a*6703494707970582.0, GMP_RNDD, "1.3853ee10c9c98@13");
434: check4(a*8010323124937260.0, GMP_RNDD, "1.556abe212b56e@13");
435: check4(a*2.0*8010776873384260.0, GMP_RNDD, "1.e2d9a51977e6d@13");
436:
437: return 0;
438: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tsqrt.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