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