Annotation of OpenXM_contrib/gmp/mpfr/tests/tagm.c, Revision 1.1.1.1
1.1 ohara 1: /* Test file for mpfr_agm.
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 "mpfr.h"
28: #include "mpfr-test.h"
29:
30: double drand_agm _PROTO((void));
31: double dagm _PROTO((double, double));
32: void check4 _PROTO((double, double, mp_rnd_t, double));
33: void check_large _PROTO((void));
34: void slave _PROTO((int, int));
35:
36: double
37: drand_agm (void)
38: {
39: double d; long int *i;
40:
41: i = (long int*) &d;
42: do {
43: i[0] = LONG_RAND();
44: i[1] = LONG_RAND();
45: } while ((d<1e-153)||(d>1e153)); /* to avoid underflow or overflow
46: in double calculus in sqrt(u*v) */
47:
48: return d;
49: }
50:
51: double
52: dagm (double a, double b)
53: {
54: double u, v, tmpu, tmpv;
55:
56: if ((isnan(a))||(isnan(b)))
57: return a+b;
58:
59: tmpv=MAX(a,b);
60: tmpu=MIN(a,b);
61:
62: do
63: {
64: u=tmpu;
65: v=tmpv;
66: tmpu=sqrt(u*v);
67: tmpv=(u+v)/2.0;
68: }
69: while (!(((tmpu==u)&&(tmpv==v))||(ulp(u,v)==0)));
70:
71: /* printf("difference : %i ulp\n",ulp(u,v)); */
72: return u;
73: }
74:
75: #define check(a,b,r) check4(a,b,r,0.0)
76:
77: void
78: check4 (double a, double b, mp_rnd_t rnd_mode, double res1)
79: {
80: mpfr_t ta, tb, tres;
81: double res2;
82: int ck=0;
83:
84: mpfr_init2(ta, 53);
85: mpfr_init2(tb, 53);
86: mpfr_init2(tres, 53);
87:
88: mpfr_set_d(ta, a, rnd_mode);
89: mpfr_set_d(tb, b, rnd_mode);
90:
91: mpfr_agm(tres, ta, tb, rnd_mode);
92: #ifdef MPFR_HAVE_FESETROUND
93: mpfr_set_machine_rnd_mode(rnd_mode);
94: #endif
95:
96: if (res1==0.0) res1=dagm(a,b); else ck=1;
97: if (ck==0) printf("%1.20e\n", res1);
98: res2 = mpfr_get_d1 (tres);
99:
100: if (ck && res1!=res2 && (!isnan(res1) || !isnan(res2))) {
101: printf("mpfr_agm failed for a=%1.20e, b=%1.20e, rnd_mode=%d\n",a,b,rnd_mode);
102: printf("expected result is %1.20e, got %1.20e (%d ulp)\n",res1,res2,
103: ulp(res2,res1));
104: /*exit(1);*/
105: }
106: mpfr_clear(ta); mpfr_clear(tb); mpfr_clear(tres);
107: }
108:
109: void
110: check_large (void)
111: {
112: mpfr_t a, b, agm;
113:
114: mpfr_init2(a, 82); mpfr_init2(b, 82); mpfr_init2(agm, 82);
115: mpfr_set_ui(a, 1, GMP_RNDN);
116: mpfr_set_str_raw(b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39");
117: mpfr_agm(agm, a, b, GMP_RNDN);
118: mpfr_set_str_raw(a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4");
119: if (mpfr_cmp(agm, a)) {
120: fprintf(stderr, "mpfr_agm failed for precision 82\n"); exit(1);
121: }
122: mpfr_clear(a); mpfr_clear(b); mpfr_clear(agm);
123: }
124:
125: void
126: slave (int N, int p)
127: {
128: int i;
129: double a,b;
130: mpfr_t ta, tb, tres;
131:
132: mpfr_init2(ta, 53);
133: mpfr_init2(tb, 53);
134: mpfr_init2(tres, p);
135: for(i=0;i<N;i++) {
136: a = drand_agm();
137: b = drand_agm();
138: mpfr_set_d(ta, a, GMP_RNDN);
139: mpfr_set_d(tb, b, GMP_RNDN);
140: mpfr_agm(tres, ta, tb, LONG_RAND() % 4 );
141: }
142: mpfr_clear(ta); mpfr_clear(ta); mpfr_clear(tres);
143: printf("fin\n");
144: }
145:
146:
147: int
148: main (int argc, char* argv[])
149: {
150: int N;
151:
152: SEED_RAND (time(NULL));
153:
154: if (argc==3) { /* tagm N p : N calculus with precision p*/
155: printf("Doing %d random tests in %d precision\n",atoi(argv[1]),atoi(argv[2]));
156: slave(atoi(argv[1]),atoi(argv[2]));
157: return 0;
158: }
159:
160: if (argc==2) { /* tagm N: N tests with random double's */
161: int i;
162: double a,b;
163:
164: N = atoi(argv[1]);
165: for (i=0;i<N;i++) {
166: a = drand();
167: b = drand();
168: check(a, b, LONG_RAND() % 4);
169: }
170: return 0;
171: }
172: else {
173: check_large();
174: check4(2.0, 1.0, GMP_RNDN, 1.45679103104690677029);
175: check4(6.0, 4.0, GMP_RNDN, 4.94936087247260925182);
176: check4(62.0, 61.0, GMP_RNDN, 6.14989837188450749750e+01);
177: check4(0.5, 1.0, GMP_RNDN, 7.28395515523453385143e-01);
178: check4(1.0, 2.0, GMP_RNDN, 1.45679103104690677029);
179: check4(234375765.0, 234375000.0, GMP_RNDN, 2.3437538249984395504e8);
180: check4(8.0, 1.0, GMP_RNDU, 3.615756177597362786);
181: check4(1.0, 44.0, GMP_RNDU, 1.33658354512981247808e1);
182: check4(1.0, 3.7252902984619140625e-9, GMP_RNDU, 7.55393356971199025907e-02);
183: }
184:
185: /* TODO : tests des infinis dans tagm.c */
186:
187: return 0;
188: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tagm.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