1 /* -*-C-*- tpower.c */
2
3 #include "elefunt.h"
4
5 /*# program to test power function (**)
6 #
7 # data required
8 #
9 # none
10 #
11 # subprograms required from this package
12 #
13 # machar - an environmental inquiry program providing
14 # information on the floating-point arithmetic
15 # system. note that the call to machar can
16 # be deleted provided the following six
17 # parameters are assigned the values indicated
18 #
19 # ibeta - the radix of the floating-point system
20 # it - the number of base-ibeta digits in the
21 # significand of a floating-point number
22 # minexp - the largest in magnitude negative
23 # integer such that float(ibeta)**minexp
24 # is a positive floating-point number
25 # maxexp - the largest positive integer exponent
26 # for a finite floating-point number
27 # xmin - the smallest non-vanishing floating-point
28 # power of the radix
29 # xmax - the largest finite floating-point
30 # number
31 #
32 # ran(k) - a function subprogram returning random real
33 # numbers uniformly distributed over (0,1)
34 #
35 #
36 # standard fortran subprograms required
37 #
38 # abs, alog, amax1, exp, float, sqrt
39 #
40 #
41 # latest revision - december 6, 1979
42 #
43 # author - w. j. cody
44 # argonne national laboratory
45 #
46 #**********************************************************************/
47
48 void
tpower(VOID_ARG)49 tpower(VOID_ARG)
50 {
51 int i,
52 ibeta,
53 iexp,
54 irnd,
55 it,
56 j,
57 k1,
58 k2,
59 k3,
60 machep,
61 maxexp,
62 minexp,
63 n,
64 negep,
65 ngrd;
66
67 dp_t
68 eps,
69 epsneg,
70 xmax,
71 xmin;
72
73 volatile dp_t
74 a,
75 ait,
76 albeta,
77 alxmax,
78 b,
79 beta,
80 c,
81 del,
82 dely,
83 onep5,
84 r6,
85 r7,
86 scale,
87 xl,
88 xn,
89 xsq,
90 x1,
91 y,
92 y1,
93 y2,
94 z,
95 zz;
96
97 dp_t
98 w,
99 x;
100
101 /*******************************************************************/
102
103 (void)ranset(initseed());
104 machar(&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
105 &maxexp, &eps, &epsneg, &xmin, &xmax);
106 beta = TO_FP_T(ibeta);
107 albeta = ALOG(beta);
108 ait = TO_FP_T(it);
109 alxmax = ALOG(xmax);
110 onep5 = (TWO + ONE) / TWO;
111 scale = ONE;
112 j = (it + 1) / 2;
113 for (i = 1; i <= j; i++) /* do i = 1, j */
114 scale = scale * beta;
115 a = ONE / beta;
116 b = ONE;
117 c = -AMAX1(alxmax, -ALOG(xmin)) / ALOG(100e0);
118 dely = -c - c;
119 n = maxtest();
120 xn = TO_FP_T(n);
121 y = ZERO;
122 y1 = ZERO;
123 zz = ZERO;
124
125 /* random argument accuracy tests */
126
127 for (j = 1; j <= 4; j++)
128 { /* do j = 1, 4 */
129 k1 = 0;
130 k3 = 0;
131 x1 = ZERO;
132 r6 = ZERO;
133 r7 = ZERO;
134 del = (b - a) / xn;
135 xl = a;
136
137 for (i = 1; i <= n; i++)
138 { /* DO I = 1, N */
139 x = del * RAN() + xl;
140 if (j == 1)
141 {
142 zz = POW(x, ONE); /* zz=x**ONE */
143 z = x;
144 }
145 else
146 {
147 w *= scale;
148 w = STORE(&w);
149 x += w;
150 x = STORE(&x);
151 x -= w;
152 x = STORE(&x);
153 xsq = x * x;
154 if (j != 4)
155 {
156 zz = POW(xsq, onep5); /* zz = xsq ** onep5; */
157 z = x * xsq;
158 }
159 else
160 {
161 y = dely * RAN() + c;
162 y2 = (y / TWO + y) - y;
163 y = y2 + y2;
164 z = POW(x, y); /* z=x ** y; */
165 zz = POW(xsq, y2); /* zz=xsq ** y2; */
166 }
167 }
168 w = ONE;
169 if (z != ZERO) /* */
170 w = (z - zz) / z;
171 if (w > ZERO)
172 k1 = k1 + 1;
173 if (w <= ZERO)
174 k3 = k3 + 1;
175 w = ABS(w);
176 if (w > r6)
177 {
178 r6 = w;
179 x1 = x;
180 if (j == 4)
181 y1 = y;
182 }
183 r7 = r7 + w * w;
184 xl = a + TO_FP_T(i) * del; /* OLD: xl = xl + del (bad for large n) */
185 }
186
187 k2 = n - k3 - k1;
188 r7 = SQRT(r7 / xn);
189 if (j == 1)
190 {
191 (void)printf("1TEST OF X**1.0 VS X\n\n\n");
192 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
193 (void)printf(" (%15.4e,%15.4e)\n\n\n", a, b);
194 (void)printf(" X**1.0 WAS LARGER%6d TIMES,\n", k1);
195 (void)printf(" AGREED%6d TIMES, AND\n", k2);
196 (void)printf(" WAS SMALLER%6d TIMES.\n\n\n", k3);
197 }
198 else if (j != 4)
199 {
200 (void)printf("1TEST OF XSQ**1.5 VS XSQ*X\n\n\n");
201 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
202 (void)printf(" (%15.4e,%15.4e)\n\n\n", a, b);
203 (void)printf(" X**1.5 WAS LARGER%6d TIMES,\n", k1);
204 (void)printf(" AGREED%6d TIMES, AND\n", k2);
205 (void)printf(" WAS SMALLER%6d TIMES.\n\n\n", k3);
206 }
207 else
208 {
209 (void)printf("1TEST OF X**Y VS XSQ**(Y/2)\n\n\n");
210 w = c + dely;
211 (void)printf(" %6d RANDOM ARGUMENTS WERE TESTED FROM THE REGION\n", n);
212 (void)printf(" X IN (%15.4e,%15.4e), Y IN (%15.4e,%15.4e)\n\n\n",
213 a, b, c, w);
214 (void)printf(" X**Y WAS LARGER%6d TIMES,\n", k1);
215 (void)printf(" AGREED%6d TIMES, AND\n", k2);
216 (void)printf(" WAS SMALLER%6d TIMES.\n\n\n", k3);
217 }
218 (void)printf(" THERE ARE %4d BASE %4d SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER\n\n\n",
219 it, ibeta);
220 w = -999.0e0;
221 if (r6 != ZERO)
222 w = ALOG(ABS(r6)) / albeta;
223 if (j != 4)
224 {
225 (void)printf(" THE MAXIMUM RELATIVE ERROR OF %15.4e = %4d ** %7.2f\n",
226 r6, ibeta, w);
227 (void)printf(" OCCURRED FOR X =%17.6e\n", x1);
228 }
229 if (j == 4)
230 {
231 (void)printf(" THE MAXIMUM RELATIVE ERROR OF %15.4e = %4d ** %7.2f\n",
232 r6, ibeta, w);
233 (void)printf(" OCCURRED FOR X =%17.6e Y =%17.6e\n", x1, y1);
234 }
235 w = AMAX1(ait + w, ZERO);
236 (void)printf(
237 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2f\n\n\n",
238 ibeta, w);
239 w = -999.0e0;
240 if (r7 != ZERO)
241 w = ALOG(ABS(r7)) / albeta;
242 (void)printf(
243 " THE ROOT MEAN SQUARE RELATIVE ERROR WAS %15.4e = %4d **%7.2f\n",
244 r7, ibeta, w);
245 w = AMAX1(ait + w, ZERO);
246 (void)printf(
247 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2f\n\n\n",
248 ibeta, w);
249 if (j != 1)
250 {
251 b = 10.0e0;
252 a = 0.01e0;
253 if (j != 3)
254 {
255 a = ONE;
256 b = EXP(alxmax / 3.0e0);
257 }
258 }
259 }
260
261 /* special tests */
262
263 (void)printf("1SPECIAL TESTS\n\n\n");
264 (void)printf(" THE IDENTITY X ** Y = (1/X) ** (-Y) WILL BE TESTED.\n\n");
265 (void)printf(" X Y (X**Y-(1/X)**(-Y)) / X**Y\n\n");
266 b = 10.0e0;
267
268 for (i = 1; i <= 5; i++)
269 { /* DO I = 1, 5 */
270 x = RAN() * b + ONE;
271 y = RAN() * b + ONE;
272 z = POW(x, y); /* z=x ** y */
273 zz = POW((ONE / x), -y);/* zz=(ONE/x) ** (-y) */
274 w = (z - zz) / z;
275 (void)printf("%15.7e%15.7e %15.7e\n\n", x, y, w);
276 }
277
278 /* test of error returns */
279
280 (void)printf("1TEST OF ERROR RETURNS\n\n\n");
281 x = beta;
282 y = TO_FP_T(minexp);
283 (void)printf(" (%14.7e) ** (%14.7e) WILL BE COMPUTED.\n", x, y);
284 (void)printf(" THIS SHOULD NOT TRIGGER AN ERROR MESSAGE\n\n\n");
285 fflush(stdout);
286 errno = 0;
287 z = POW(x, y); /* Z = X ** Y */
288 if (errno)
289 perror("POW()");
290 (void)printf(" THE VALUE RETURNED IS %15.4e\n\n\n\n", z);
291 y = TO_FP_T(maxexp - 1);
292 (void)printf(" (%14.7e) ** (%14.7e) WILL BE COMPUTED.\n", x, y);
293 (void)printf(" THIS SHOULD NOT TRIGGER AN ERROR MESSAGE\n\n\n");
294 fflush(stdout);
295 errno = 0;
296 z = POW(x, y); /* z = X ** Y */
297 if (errno)
298 perror("POW()");
299 (void)printf(" THE VALUE RETURNED IS %15.4e\n\n\n\n", z);
300 x = ZERO;
301 y = TWO;
302 (void)printf(" (%14.7e) ** (%14.7e) WILL BE COMPUTED.\n", x, y);
303 (void)printf(" THIS SHOULD NOT TRIGGER AN ERROR MESSAGE\n\n\n");
304 fflush(stdout);
305 errno = 0;
306 z = POW(x, y); /* z = X ** Y */
307 if (errno)
308 perror("POW()");
309 (void)printf(" THE VALUE RETURNED IS %15.4e\n\n\n\n", z);
310 x = -y;
311 y = ZERO;
312 (void)printf(" (%14.7e) ** (%14.7e) WILL BE COMPUTED.\n", x, y);
313 (void)printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
314 fflush(stdout);
315 errno = 0;
316 z = POW(x, y); /* z = X ** Y */
317 if (errno)
318 perror("POW()");
319 (void)printf(" THE VALUE RETURNED IS %15.4e\n\n\n\n", z);
320 y = TWO;
321 (void)printf(" (%14.7e) ** (%14.7e) WILL BE COMPUTED.\n", x, y);
322 (void)printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
323 fflush(stdout);
324 errno = 0;
325 z = POW(x, y); /* z = X ** Y */
326 if (errno)
327 perror("POW()");
328 (void)printf(" THE VALUE RETURNED IS %15.4e\n\n\n\n", z);
329 x = ZERO;
330 y = ZERO;
331 (void)printf(" (%14.7e) ** (%14.7e) WILL BE COMPUTED.\n", x, y);
332 (void)printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
333 fflush(stdout);
334 errno = 0;
335 z = POW(x, y); /* z = X ** Y */
336 if (errno)
337 perror("POW()");
338 (void)printf(" THE VALUE RETURNED IS %15.4e\n\n\n\n", z);
339 (void)printf(" THIS CONCLUDES THE TESTS\n");
340 }
341