1 /* -*-C-*- talog.c */
2
3 #include "elefunt.h"
4
5 /*# program to test alog
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 four
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 # xmin - the smallest non-vanishing floating-point
23 # power of the radix
24 # xmax - the largest finite floating-point no.
25 #
26 # ran(k) - a function subprogram returning random real
27 # numbers uniformly distributed over (0,1)
28 #
29 #
30 # standard fortran subprograms required
31 #
32 # abs, alog, alog10, amax1, float, sign, sqrt
33 #
34 #
35 # latest revision - december 6, 1979
36 #
37 # author - w. j. cody
38 # argonne national laboratory
39 #
40 #*/
41
42 void
talog(VOID_ARG)43 talog(VOID_ARG)
44 {
45 int i,
46 ibeta,
47 iexp,
48 irnd,
49 it,
50 j,
51 k1,
52 k2,
53 k3,
54 machep,
55 maxexp,
56 minexp,
57 n,
58 negep,
59 ngrd;
60
61 qp_t
62 eps,
63 epsneg,
64 xmax,
65 xmin;
66
67 volatile qp_t
68 a,
69 ait,
70 albeta,
71 b,
72 beta,
73 c,
74 del,
75 eight,
76 half,
77 r6,
78 r7,
79 tenth,
80 w,
81 x,
82 xl,
83 xn,
84 x1,
85 y,
86 z,
87 zz;
88
89 /*******************************************************************/
90
91 (void)qranset(initseed());
92 macharl(&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
93 &maxexp, &eps, &epsneg, &xmin, &xmax);
94
95 beta = TO_FP_T(ibeta);
96 albeta = ALOG(beta);
97 ait = TO_FP_T(it);
98 j = it / 3;
99 half = 0.5e+00L;
100 eight = 8.0e+00L;
101 tenth = 0.1e+00L;
102 c = ONE;
103
104 for (i = 1; i <= j; ++i)
105 c = c / beta;
106
107 b = ONE + c;
108 a = ONE - c;
109 n = maxtest();
110 xn = TO_FP_T(n);
111
112 /* random argument accuracy tests */
113
114 for (j = 1; j <= 4; ++j)
115 {
116 k1 = 0;
117 k3 = 0;
118 x1 = ZERO;
119 r6 = ZERO;
120 r7 = ZERO;
121 del = (b - a) / xn;
122 xl = a;
123
124 for (i = 1; i <= n; ++i)
125 {
126 x = del * RAN() + xl;
127 if (j == 1)
128 {
129 y = (x - half);
130 y -= half;
131 zz = ALOG(x);
132 z = ONE / 3.0e+00L;
133 z = y * (z - y / 4.0e+00L);
134 z = (z - half) * y * y + y;
135 }
136 else if (j == 2)
137 {
138 x += eight;
139 x -= eight;
140 y = x + x / 16.0e+00L;
141 z = ALOG(x);
142 zz = ALOG(y) - 7.77468164348425806061320404202632862024751447237708145e-05L;
143 zz = zz - 31.0e+00L / 512.0e+00L;
144 }
145 else
146 if (j != 3)
147 {
148 z = ALOG(x * x);
149 zz = ALOG(x);
150 zz = zz + zz;
151 }
152 else
153 {
154 x += eight;
155 x -= eight;
156 y = x + x * tenth;
157 z = ALOG10(x);
158 zz = ALOG10(y) - 3.770601582250407501999712430242417067021904664530945966e-04L;
159 zz = zz - 21.0e+00L / 512.0e+00L;
160 }
161 w = ONE;
162 if (z != ZERO)
163 w = (z - zz) / z;
164 z = SIGN(w, z);
165 if (z > ZERO)
166 k1 = k1 + 1;
167 if (z < ZERO)
168 k3 = k3 + 1;
169 w = ABS(w);
170 if (w > r6)
171 {
172 r6 = w;
173 x1 = x;
174 }
175 r7 = r7 + w * w;
176 xl = a + TO_FP_T(i) * del; /* OLD: xl = xl + del (bad for large n) */
177 }
178
179 k2 = n - k3 - k1;
180 r7 = SQRT(r7 / xn);
181 if (j == 1)
182 (void)printf("1TEST OF ALOG(X) VS T.S. EXPANSION OF ALOG(1+Y)\n\n\n");
183 else if (j == 2)
184 (void)printf("1TEST OF ALOG(X) VS ALOG(17X/16)-ALOG(17/16)\n\n\n");
185 else if (j == 3)
186 (void)printf("1TEST OF ALOG10(X) VS ALOG10(11X/10)-ALOG10(11/10)\n\n\n");
187 else if (j == 4)
188 (void)printf("1TEST OF ALOG(X*X) VS 2 * LOG(X)\n\n\n");
189 if (j == 1)
190 {
191 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n",n);
192 (void)printf(" (1-EPS,1+EPS), WHERE EPS =%15.4Le\n\n\n", c);
193 }
194 else
195 {
196 (void)printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
197 (void)printf(" (%15.4Le,%15.4Le)\n\n\n", a, b);
198 }
199 if (j != 3)
200 {
201 (void)printf(" ALOG(X) WAS LARGER%6d TIMES,\n", k1);
202 (void)printf(" AGREED%6d TIMES, AND\n", k2);
203 (void)printf(" WAS SMALLER%6d TIMES.\n\n\n", k3);
204 }
205 else
206 {
207 (void)printf(" ALOG10(X) WAS LARGER%6d TIMES,\n", k1);
208 (void)printf(" AGREED%6d TIMES, AND\n", k2);
209 (void)printf(" WAS SMALLER%6d TIMES.\n\n\n", k3);
210 }
211 (void)printf(" THERE ARE%4d BASE%4d SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER\n\n\n",
212 it, ibeta);
213 w = -999.0e+00L;
214 if (r6 != ZERO)
215 w = ALOG(ABS(r6)) / albeta;
216 (void)printf(" THE MAXIMUM RELATIVE ERROR OF%15.4Le = %4d **%7.2Lf\n",
217 r6, ibeta, w);
218 (void)printf(" OCCURRED FOR X =%17.6Le\n", x1);
219 w = AMAX1(ait + w, ZERO);
220 (void)printf(
221 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2Lf\n\n\n",
222 ibeta, w);
223 w = -999.0e+00L;
224 if (r7 != ZERO)
225 w = ALOG(ABS(r7)) / albeta;
226 (void)printf(" THE ROOT MEAN SQUARE RELATIVE ERROR WAS%15.4Le = %4d **%7.2Lf\n",
227 r7, ibeta, w);
228 w = AMAX1(ait + w, ZERO);
229 (void)printf(
230 " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS%7.2Lf\n\n\n",
231 ibeta, w);
232 if (j <= 1)
233 {
234 a = SQRT(half);
235 b = 15.0e+00L / 16.0e+00L;
236 }
237 else if (j > 2)
238 {
239 a = 16.0e+00L;
240 b = 240.0e+00L;
241 }
242 else
243 {
244 a = SQRT(tenth);
245 b = 0.9e+00L;
246 }
247 }
248
249 /* special tests */
250
251 (void)printf("1SPECIAL TESTS\n\n\n");
252 (void)printf(" THE IDENTITY ALOG(X) = -ALOG(1/X) WILL BE TESTED.\n\n");
253 (void)printf(" X F(X) + F(1/X)\n\n");
254
255 for (i = 1; i <= 5; ++i)
256 {
257 x = RAN();
258 x = x + x + 15.0e+00L;
259 y = ONE / x;
260 z = ALOG(x) + ALOG(y);
261 (void)printf("%15.7Le%15.7Le\n\n", x, z);
262 }
263
264 (void)printf("\n\n TEST OF SPECIAL ARGUMENTS\n\n\n");
265 x = ONE;
266 y = ALOG(x);
267 (void)printf(" ALOG(1.0) = %15.7Le\n\n\n", y);
268 x = xmin;
269 y = ALOG(x);
270 (void)printf(" ALOG(XMIN) = ALOG(%15.7Le) = %15.7Le\n\n\n", x, y);
271 x = xmax;
272 y = ALOG(x);
273 (void)printf(" ALOG(XMAX) = ALOG(%15.7Le) = %15.7Le\n\n\n", x, y);
274
275 /* test of error returns */
276
277 (void)printf("1TEST OF ERROR RETURNS\n\n\n");
278
279 x = -2.0e+00L;
280 (void)printf(" ALOG WILL BE CALLED WITH THE ARGUMENT%15.4Le\n",x);
281 (void)printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
282 fflush(stdout);
283 errno = 0;
284 y = ALOG(x);
285 if (errno)
286 perror("ALOG()");
287 (void)printf(" ALOG RETURNED THE VALUE%15.4Le\n\n\n", y);
288
289 x = ZERO;
290 (void)printf(" ALOG WILL BE CALLED WITH THE ARGUMENT%15.4Le\n",x);
291 (void)printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
292 fflush(stdout);
293 errno = 0;
294 y = ALOG(x);
295 if (errno)
296 perror("ALOG()");
297 (void)printf(" ALOG RETURNED THE VALUE%15.4Le\n\n\n\n", y);
298
299 (void)printf(" THIS CONCLUDES THE TESTS\n");
300 }
301